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Abstract 


In this thesis, we investigate two problems of entangled branched polymers, i.e., 
the numerical solutions of the arm-retraction problem for well entangled star arms 
and the relaxation behaviours of branched polymers with different architectures. 
For the first problem, the arm retraction dynamics is studied using both the one- 
dimensional Rouse chain model and the slip-spring model by an advanced numerical 
method for the first-passage time problems, namely the forward flux sampling (FFS) 
method. In the one-dimensional Rouse chain model, we measured the first-passage 
time that the arm free end extends to a distance away from the origin, showing 
that the mean first-passage time is getting shorter if the Rouse chain is represented 
by more beads. The simulation results validate the prediction of an asymptotic 
solution for the multi-dimensional first-passage problem, which suggests the arm 
retraction time is much shorter than the prediction of the Milner-McLeish theory 
without constraint release. Then, we implement the FFS method to the slip-spring 
model and get the relaxation spectra for different arm lengths, ranging from mildly 
entangled to well-entangled star arms. We also proposed an algorithm to extract 
the dynamic observables, i.e., the end-to-end vector and stress relaxation functions, 
from the FFS simulation results. For the second problem, we conduct a series of 
molecular dynamics (MD) simulations using high performance GPU methods on the 
mildly entangled branched polymers of different architectures, including 3-arm sym- 
metric and asymmetric stars, and H-shaped polymers. The slip-spring model, whose 


parameters are carefully calibrated according to the MD results of linear chains, is 


il 


also implemented to predict the relaxation behaviours of the branched polymers. We 
present a detailed analysis on the arm end-to-end vector relaxation functions and 
the monomer mean-squared displacements. By comparing the MD and slip-spring 
model simulation results, we propose a slip-link “hopping” mechanism, which ac- 
counts for the behaviour that the entanglements can pass through the branch point 


when the third arm is disentangled. 
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Chapter 1 


Introduction 


1.1 Background 


The annual worldwide production of synthetic polymers has reached several mil- 
lion tons and still undergoes continuous growth. Before being shaped into com- 
mercial plastic or rubber products, the raw materials are processed in molten or 
liquid state, thus our understanding of their rheology is crucial to enhancing their 
processing properties. Driven by industrial interests and scientific curiosity, both 
experimental and theoretical research on the rheological properties of bulk polymer 
fluids has experienced a fast development in the past half century. It is known that 
polymeric liquids have very complicated rheological and thus processing properties 
due to the hierarchical relaxation behaviours at different time and length scales. 
These behaviours are governed by polymer architectures and distribution of molec- 
ular weights. A typical example is that the presence of a small amount of long 
chain branching (LCB) in commercial polymers can dramatically affect rheolog- 
ical behaviour, e.g., giving them a much higher zero-shear viscosity 7 than in linear 
systems with the same molecular weights. Understanding the relationship between 
LCB and rheology could simultaneously solve both the direct and inverse problems: 


to predict the rheological and thus the processing properties of a given polymer 
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system and using rheological measurements to deduce molecular weight distribu- 
tion and identify LCB. Such relationship, however, remains to be one of the most 
challenging problems in polymer physics. 

A typical polymer chain is constituted by repeated chemical units, called monomers, 
connected via covalent bonds. When polymers are in bulk, the polymer melts ex- 
hibit a fascinating property called “viscoelasticity”. The stress of the viscoelastic 
materials neither simply depends on the strain nor the strain rate, but is a function 
of the deformation history [I] [2]. For example, in response to a small step strain y, 
the stress o of of an isotropic polymer liquid first raises to a level proportional to the 
strain, then slowly decays over a long period. Such behaviour showing both elastic- 
ity and viscosity could be qualitatively described by the general Maxwell model, 
which is constructed by connecting a set of linear springs and dashpots in series [3]. 


The stress modulus of the this mechanical model, G(t), is a sum of exponentials: 
G(t) = S Gi exp (—t/7;) 
i=1 


where G; and 7; are the moduli and the relaxation times of the Maxwell modes. 
The stress relaxation modulus exhibits an exponential decay beyond the longest 
relaxation time 7,,, which is also referred to as the terminal relaxation time Tq in 
literature due to its significance. 

The complexity of polymer dynamics is far beyond the description of the phe- 
nomenological mechanical models, which can be reflected in the rheological prop- 
erties of polymers [4H7], and thus other measurements analogous to rheology, such 
as dielectric [8], optical [9] and diffusional properties. A typical example 
is the dependence of rheological properties on the molecular weight, MM. For linear 
polymer melts, the viscosity 7) increases linearly with M at low molecular weights, 


but exhibits power law growth of M?* when M exceeds a critical value, M: 


M, M<™M, 


No ~ 
Mee. MSM. 


CHAPTER 1. INTRODUCTION 3 


log(G(t)) 


Tq 


e log(t) 


Figure 1.1: The typical relaxation modulus G(t) of a well-entangled linear polymer 


melt after a step-strain. 


Fig. sketches the stress relaxation modulus of a linear polymer melt with 
M > M,. At short time scales, the stress relaxation moduli cannot be distinguished 
between different molecular weights. After a certain time 7, G(t) for polymers with 
a large molecular weight barely decays, exhibiting a plateau with a M-independent 
modulus, G,. Such plateau regime cannot be observed in low-M polymer systems, 
but is analogous to the plateau modulus of a cross-linked polymer network with a 
strand molecular weight I, between the cross-links [13]. The critical molecular 
weight MM, is approximately 2M, for all amorphous melts independent of their chem- 
istry. After the plateau regime, G(t) decays exponentially with the characteristic 
time Tg which grows with M as M?*. It is well known that the remarkable differ- 
ence in the rheological behaviours of linear melts with different molecular weights is 
attributed to “entanglements”, which are essentially the topological interactions 
between neighboring polymers. In a melt of short-chain polymers, a probe chain 
diffuses freely as in a frictional fluid, whose dynamics could be well interpreted by 
the Rouse model [14]. In the long-chain polymer melts, the probe polymer chain 
feels the topological constraint, by which the chains cannot cross each other, thus 
lateral diffusion against such constraints is prohibited, leading to a drastic slowing 


down on diffusion. 
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The most successful theoretical framework for entangled polymers is the famous 
tube model originally proposed by de Gennes, Doi and Edwards (T5Hi7]_, 
based on which theories for entangled dynamics have developed for half century, 
accounting for extensive experimental data. In the fundamental assumption of the 
tube model, the topological constraints imposed by surrounding chains act to confine 
the probe chain in a tube-like region. The contour line of the tube is called the prim- 
itive path, upon which the projected chain is called the primitive chain. The 
monomer diffusion perpendicular to the primitive path is restricted in the tube-like 
region with a diameter a. The chain can only diffuse curvilinearly along the primitive 
path at time scales longer than the entanglement time 7, which is the time taken by 
an entanglement segment to diffuse a distance of its own size a, i.e., the tube diam- 
eter. This type of diffusion is known as “reptation” [16]. “Reptation”, which means 
“creeping” in Latin, represents the snake-like motion of a polymer chain in its own 
tube. On the length scales larger than the tube diameter a, reptation is equivalent 
to the one-dimensional diffusion of the primitive chain: the central sections of the 
chain follow the tube contour while the chain ends explore in the melt to create new 
tube segments. With very limited parameters, the original tube model successfully 
captures the qualitative features of entanglement dynamics of linear chains, but fails 
to provide quantitative predictions. For example, the predicted dependence of vis- 
cosity on molecular weight is 7) ~ M*° rather than the experimentally observed 
M?*. Such discrepancies arise from the oversimplified assumptions made in the 
original tube model, such as the inextensibility of the primitive chain. 

In order to provide quantitative predictions, the tube model needs to be im- 
proved by incorporating additional relaxation modes, such as the contour length 
fluctuation (CLF) due to the fluctuations of tube length around the average value, 
and the constraint release (CR) due to the exchange of surrounding chains. In 
the picture of reptation plus CLF, Doi deduced the 3.4 power law as well as 


the particular range it accounts for, which quantitatively predicted the experimen- 
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tal observation that the power law for very long chains would deviate from 3.4 to 
3.0. The CR mechanism, which is extraordinarily important in polydisperse systems 
where the molecular weights of different species are significantly separated from each 
other, represents the effect brought by the exchange of the surrounding chains on 
the tube confinement. In a CR event, a probe chain is permitted to move transverse 
to the confining tube when some neighbouring chains move away by reptation or 
CLF. Graessley |6} assumes that the exchange of the surroundings chains does 
not modify the tube diameter, but results in a rearrangement of tube segments 
of the probe chain. So the tube becomes a random walk and thus amounts to a 
Rouse-like motion, which is also referred to as the constraint-release Rouse (CR 
Rouse) motion. Cloizeaux proposed a “double-reptation” approximation and 
inferred the stress relaxation modulus Gt) as a quadratic function of the tube sur- 
vival fraction p(t). This “double-reptation” picture works best for the polydisperse 
systems with broad molecular weight distributions, but show limitation in binary 
blends of chains with greatly separated lengths. Rubinstein [2123] considered the 
broad distribution of the CR rates, and proposed a self-consistent theory for binary 
blends, where the stress relaxation modulus G(t) is in a product form of the tube 
survival fraction s(t) and a Rouse-like relaxation function R(t) of the tube repre- 
senting the CR process. By incorporating CLF and CR, the tube model explained 
the additional stress relaxed at the reptation time of short chains in comparison 
with the pure reptation model. A quantitative tube theory for well-entangled linear 
chain systems was proposed by Likhtman and McLeish [24], who combined the main 
relaxation mechanisms in a self-consistent manner. 

As anticipated, a simple change in the polymer architectures, such as adding a 
branch to a linear chain to create a 3-arm star polymer could completely change 
the dynamics. In a f-arm star polymer (f > 3), reptation is highly suppressed 
because entropically it is extremely unfavourable to drag f — 1 arms into one tube. 


Thus the branch point would be localized in a volume of tube diameter size before 
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the arms are relaxed by CLF and CR. In star polymers, the CLF mechanism is 
also termed as “arm retraction”, in which the arm repeatedly retracts inwards 
along the primitive path toward the branch point to release original tube segment 
and stretches out to create new tube segments. de Gennes indicated that arm 
retraction in a fixed network is an exponentially-rare event, whose probability was 
later calculated by Kuzzu and Doi |26]. Based on their works, Person and Helfand 
presented a refined theory for arm retraction, which assumes the entropically 
unfavoured retraction is a thermally activated process in a quadratic potential field, 
U(s) = vkgT(M/M.,)s”, where s represents the relaxed fraction of the primitive 
path changing from zero to unity during the retraction, vy & 0.6 is a numerical 
factor obtained by fitting experimental results. The Person-Helfand theory confirms 
the experimental observation that the viscosity of star polymers grows exponentially 
with the increasing molecular weight of the arm, M,. However, the numerical factor 
v = 0.6 is in contrast to the CLF model which provides v = 15/8 [17]. 

Ball and McLeish [28] realized that the discrepancy in the value of v can be 
attributed to the absence of CR mechanism in the Person-Helfand theory. By in- 
voking a “tube-dilation” hypothesis [29], they proposed a well-known concept called 
“dynamic tube dilation” (DTD), which incorporates the CR mode into the tube 
model for branched polymers. The key assumption of DTD is that the relaxed arm 
segments could be treated as solvent for unrelaxed materials, such that the arm 
would retract in a gradually widening tube. The molecular weight of the dilated 
tube segment subjects to a dilation exponent a: M,.(s) = Moo(1 — s)~*%, where 
M.o is the molecular weight of the undilated tube segment and a is either 4/3 or 
1 whose exact value should depend on the physical nature of entanglements that 
is not yet fully resolved. When taking a = 1, the effective potential Ug(s) is 
(15/8)(M,/M.)kpT (s? — 2s3/3), which gives v = 0.625 at. s = 1 in accordance with 
vy = 0.6 in Person-Helfand theory. The DTD hypothesis is expected to work better 


for branched polymers due to the well-separated relaxation timescales. In the mean- 
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time, the continuous and smooth spectrum of relaxation time validates the gradual 
dilation ansatz. Milner and McLeish developed a theory to predict the first- 
passage time of arm retraction by solving the Kramers problem of a one-bead 
linked to the branch point by a harmonic spring. The solution gives the prefactor 


of the exponential function of the relaxation time r(s): 


aie (7) Yea) 


where a scaling of (M,/M.)*”? is found. Milner-McLeish theory also considers the 
early time fluctuations where the energy barrier is smaller than kg7’, and formulates 
T(s) into a crossover function of early to late timescales. The predictions of Milner- 
McLeish theory on the loss modulus modulus G’ (w) of symmetric star polymer melts 
are in good agreement with the experimental data for the whole range of frequency 
[32] [33]. Recently, Cao et al. [34] presented an analytical solution for arm retraction 
by solving it as a multi-dimensional first-passage time problem. By including all 
Rouse modes rather than only considering the slowest one as in Milner-McLeish 
theory, the relaxation time in the absence of CR is reduced by a factor of 2/M,, 
implying a smaller prefactor for T(s). 

After the remarkable success on monodisperse symmetric stars, the Milner- 
McLeish model was then extended to other systems, such as star-linear blends 
and asymmetric stars [86]. In a blend of monodisperse linear chains and stars, the 
DTD picture fails when the linear chains are fully relaxed by reptation at their ter- 
minal time 7g, such that the tube diameter can not increase as fast as the decrease 
of unrelaxed materials in the system. Milner et al. dealt with this problem by 
assuming three relaxation stages. In the early regime t < 7g, the standard DTD 
model is applied by treating the linear chains as two-arm stars. After Tg, the tube 
of star arms undergoes CR Rouse motion with a tube segmental relaxation time 
Ta. The CR Rouse motion ends when the tube has explored the “supertube” de- 
fined by remaining entangled arm segments, after which the DTD picture resumes. 


In the simplest asymmetric stars consisting of one short arm and two long arms, 
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the terminal relaxation time of the short arm 7, is much shorter than that of the 
long arm 7}. Each time the short arm is fully retracted, the branch point is able 
to hop by a distance of a fraction of tube diameter pa with p < 1. Thus, at the 
time scales t > 7,;, the short arm can be treated as effective frictional bead with 
a larger friction at the branch point, whose diffusion coefficient is estimated using 
Einstein relation: Dy, = p?a?/(27,) [86]. More recently, the Milner-McLeish theory 
has been generalized to model branched polymers with more complicated architec- 
tures, such as H-shaped, comb, Cayley tree polymers and their general mixtures 
44|. Relaxation of these polymers proceeds hierarchically, starting from the re- 
traction of the outermost branch arms and proceeding to inner layers until the core 
of the molecules. Under this frame, several computational models were developed 
for predicting viscoelastic properties of general mixtures of branched polymers with 
different architectures, including the hierarchical model by Larson et al. [44] [45], 
the “bob” (“branch-on-branch”) model by Das et al. and the time marching 
algorithm by van Ruymbeke et al. [46H48]. These models differ in certain compu- 
tational algorithms and relaxation mechanisms, but share similar theoretical frame- 
work that incorporates reptation, contour length fluctuation (or arm retraction), 
and constraint release by dynamic tube dilation or constraint-release Rouse motion. 
A detailed comparison between hierarchical model and “bob” model was carried out 
in Ref. [44]. 

Although tube-based theories have been shown to provide semi-quantitative pre- 
dictions for rheological behaviours of various branched polymers, they are far from 
first-principle due to many mathematical and physical approximations involved. For 
example, the empirically fitting parameter p? which determines the fractional hop- 
ping distance of the branch point is found with a broad range from 1 to 1/60 for 
different asymmetric stars. For H-polymers, however, the range of p” is relatively 
narrow, roughly from 1/12 to 1/15 [36]. Whether the value of p? should be univer- 


sal or system-dependent remains unknown. In star-linear blends, the DTD picture 
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would fail when the dilation speed is higher than the diffusion rate which the chain 
segments explore the tube. Before the DTD resumes, the confining tube of arms un- 
dergoes CR Rouse motion in a “supertube” [85]. In fact, arm retraction also happens 
in this intermediate regime, but how to handle it remains a question. 

These open questions invariably point to the same origin: despite the concepts 
of entanglement and tube having been widely used in theoretical models for half 
century, there is still a lack of clear microscopic pictures about them. Computer 
simulations using a generic bead-spring polymer chain model, namely the Kremer- 
Grest (KG) model [49] makes it possible to find the highly desired picture of 
entanglement and tube. One of the remarkable accomplishments was achieved by 
Everaers et al. [51]. They conducted a “primitive-path analysis” (PPA) to count 
the number of entanglements per chain and the average number of monomers in an 
entanglement strand, and thus predicted the plateau modulus comparable to ex- 
perimental data. Recently, Likhtman and Ponmurugan employed “mean-path 
analysis” [53] to trace the entanglements by constructing the “contact map” of 
chain segments, and got the lifetime distribution and the mean-square displacements 
of persistent binary contacts between neighbouring chains or entanglements. This 
algorithm has been applied by Cao and Wang to investigate the microscopic picture 
of constraint release in symmetric star polymer melts [54]. Likhtman [55] also pro- 
posed an algorithm to construct the “tube axis” for entangled polymers, whereby 
the binary and triple entanglements could be manifested through the smoothly av- 
eraged paths of the polymers. Likewise, the assumptions in tube theory that are 
necessarily empirical and speculative could be tested microscopically, such as the 
direct visualization of branch point hopping 57]. 

In the past two decades, the slip-link or slip-spring models which work as alter- 
native descriptions for entangled polymers have captured great attention. Different 
from the tube model that treats the topological constraints given by the surround- 


ing chains in a mean-field way, this class of models assume the entanglements as 
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binary interactions between neighbouring chain segments, and thus can incorporate 
more refinements, such as the discrete description of entanglements that distribute 
along the chains with a wide range of lifetimes. The “slip-link” picture was first 
proposed by Doi in his slip-link network model, where each slip-link repre- 
sents an entanglement, and the release of an entanglement occurs only when either 
ends of the two individual chains passes through the slip-link to destroy it. The 
later slip-link based models, without exception, follow the primary picture of Doi’s 
model, but differ in modeling and algorithms. Schieber et al. [59H61] developed a 
single-chain slip-spring model, in which each slip-link represents an entanglement 
whose dynamics is governed by a stochastic equation that exchanges the number of 
monomers between neighbouring entanglement strands. The rheological properties 
can thus be obtained by tracking the location of slip-links and number of monomers 
along each entangled strand in between two adjacent slip-links. Shanbhag and Lar- 
son developed a multi-chain slip-link model by simulating an ensemble of 
primitive chains. In this slip-link model, the entanglements are modeled by pair- 
ing the slip-links on different chains, whereby the constraint release is incorporated 
naturally by destroying the slip-link pairs when either of them escapes the ends of 
involved chains by primitive path fluctuation or reptation. The multi-chain slip-link 
model developed by Masubuchi et al. [64] {65] simulates the primitive chains bundled 
by slip-links that form a 3D network, whereby this model is called primitive chain 
network (PCN) model (also implemented in the NAPLES code). The motion of 
the slip-links is governed by the force balance between entanglement strands, while 
the diffusion of the monomers along each chain contour is operated according to the 
tension distribution along the chain. The slip-link based model developed by Likht- 
man is called slip-spring model due to the slip-links are connected to a set of 
virtual springs anchored in space, which provides extra fluctuation to the location of 
entanglements. The polymer chains are modeled as Rouse chains, whose primitive 


paths are defined by the slip-links sitting on the monomers. The contour length fluc- 
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tuation is incorporated by the Rouse motion of the chains confined in the slip-links, 
while the constraint release is included in a similar way as other slip-link models. 
Based on the slip-spring model and the primitive chain network model, Langeloth 
el al. [67| recently developed another slip-spring model which employs dissipative 
particle dynamics (DPD) to simulate a 3D chain network. This model is built 
on a finer level by incorporating excluded volume as non-bonded interaction, but in 
return requires much large system and expensive computational cost. 

The aim of this thesis is to study the relaxation dynamics of branched polymers 
with the simplest architectures, such as star and H-polymers, on which a complete 
set of observables predicted by tube theories have been compared with experiments. 
Efficient methodology to investigate the arm-retraction dynamics of well entangled 
branched polymers would be developed to examine existing tube-based theories, 
such as the Milner-McLeish theory on arm retraction dynamics [34]. We will 
also explore the possibility of extending the slip-spring model to branched polymers. 
Despite its success on describing linear polymers, the slip-spring model encounters 
difficulty in extending its applications to branched polymers. The problem may lie 
in the absence of certain relaxation mechanisms close to the branch point, which 
could possibly be resolved by incorporating the hypothesis in tube theory for asym- 
metric stars, i.e., the fully retracted short arm allows the reptation of other two long 
arms. Such attempt can reversely verify the tube theory for polymers with different 
architectures. 

The content of this thesis is arranged as follows. In the remaining part of this 
chapter, we will introduce the theoretical background of polymer dynamics. In the 
second chapter, we investigate two advanced numerical methods, namely forward 
flux sampling (FFS) and weighted ensemble (WE) [70] method, and 
their applications on the multi-dimensional first-passage time problem of arm retrac- 
tion of star polymers. The FFS simulation results are compared with theoretical 


predictions from both one-dimensional and multi-dimensional solutions of arm re- 
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traction dynamics. In the third chapter, a combinational method of FFS and the 
slip-spring model is implemented to investigate the dynamics of arm-retraction in 
star polymers. With a controllable precision, this method allows direct comparison 
between the slip-spring model and the tube theory for well-entangled star polymers 
of arm length up to 16 entanglements. Moreover, a study is conducted on the ex- 
traction of experimentally measurable observables from FFS simulations, such as 
the end-to-end relaxation and stress relaxation functions. We believe this work will 
not only expand the application of FFS method to polymer dynamics by reproduc- 
ing full dynamic spectrum rather than just the first-passage time, but also to many 
other scientific areas. In the fourth chapter, we present a multi-scale computer sim- 
ulation study on the relaxation dynamics of various branched polymers, including 
symmetric and asymmetric stars, and H-polymers. The slip-spring model is updated 
for branched polymers by incorporating a parameter free mechanism, whose results 
are compared with the fully flexible Kramer-Grest model [49]. The conclusions are 


given in the last chapter. 


1.2. Polymer Chain Models 


This section introduces the static properties of polymers by focusing on the single- 
chain conformation. By comparing the static properties of the ideal chain and 
the real chain, the fundamental concept of “universality” on chain conformation is 


concluded. 


1.2.1 Ideal Chains 


Consider a chain formed by n + 1 monomers that are connected in sequence via n 
bonds, as shown in Fig. If the monomers several bonds apart do not interact 
with each other, this chain is called ideal chain. For example, in an ideal chain, 


there is no interaction between monomer 7 and monomer i+ 7 (j >> 1), even though 
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Figure 1.3: (a) Freely rotating chain model; (b) Hindered rotation model. 


they might be very close in space. The simplest ideal chain model is the freely 
jointed chain model, in which the bond length is fixed while the orientation 
of bonds is random. Based on this model, more features could be included to 
resemble the polymer chains. For example, the orientational priority of the bonds 
can be introduced by fixing the bond angle @, which is called the freely rotating 
chain model (Fig. [1.3{a)). Likewise, a torsion potential U(y) could be introduced 
to represent the steric hindrance between functional groups, which is called the 
hindered rotation model (Fig. [1-3{b)). 

An important quantity to characterize the static property of a polymer chain is 
the end-to-end vector, R,, which is the vector between the end beads, Re = r,—Yro 


or the sum of the bond vectors, R, = )7j_, r;. For an ideal chain that is long enough, 
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the average end-to-end vector (R,) is zero, but the average mean-square end-to-end 
distance is non-zero: 

(RL) = Can’ (1.1) 
where I is the bond length, C,, is called Flory’s characteristic ratio that deter- 


mines the local stiffness of the chain. C,, for different models are listed as below: 


Freely Jointed Chain Model 


7 1+cosdé 


Freely Rotating Chain Model = 
1—cosé 


Coo 


Hindered Rotation Model 


Eq. [1.JJindicates a universal property of the ideal chains: the models having various 
stiffnesses show similar behaviors at length scales larger than C’,,/. Coarse-graining 
Cy monomers into one monomer, these ideal chains are equivalent to a freely jointed 


chain with N bonds: 
(ROS N, N=n/Cy, C= Oxi. 


where the length b is called Kuhn length. 


1.2.2 Entropic Elasticity 


Consider a freely jointed chain with N + 1 monomers. The N bonds connecting 
these monomers have a fixed bond length b. The conformation of the chain can be 
visualised as a 3D random walk with a constant step length b. In each direction, 
the 3D random walker could be decomposed into 1D random walk, whose mean- 
square step length is b?/3. Such 1D random walk follows an Ornstein-Uhlenbeck 


process, thus the end-to-end distance of the chain in each direction follows Gaussian 


3 3R2, 
Pal, Rea) = Vi aang? \ one 


distribution: 
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where a is the Cartesian component. Therefore, the end-to-end vector R, of an ideal 


chain also follows Gaussian distribution, which is a product of the 3 components: 


a° Se 3R? 
Pra(N, R.) = Pra(N, Ries) PigiN, Reg) Pia, Fs) = (<7) exp (-s55) : 


(1.2) 
The distribution above indicates that the free chain is keen to collapse into coil 
with a zero end-to-end vector (chain ends tend to overlap). If the chain is stretched 


with an end-to-end vector of Re, the free energy F’ increases by 
AF = F(N,R.) — F(N,0) = -—TAS = —kpT (InQ(N, R,) — Mn Q(N,0)), 


where S is the entropy, kp is the Boltzmann constant, Q(N,R.) is the number of 
the states with end-to-end vector equal to R,. Then the statistical weight of the 
conformations with the end-to-end vector equal to R, is given by 


Q(N, R.) 


Pag(N, Re) = 
al, Re) f Q(N,R.)dR. 


Then we have, 


3 R? 
InQ(N, Re) = In(P3a(N, Re)) + Const = aE + Const. 
Therefore, the increase of the free energy AF is a quadratic function: 


_ 3kpT p> 


AP = oNPR e 


(1.3) 


which is effectively equal to a harmonic spring with an elastic constant 3kg7'/N0. 
The elasticity of a ideal chain due to the change of the entropy is called entropic 


elasticity. 


1.2.3. Real Chain 


Different from the ideal chains, a real chain has long-range interactions, which means 
that two monomers in one chain can interact even if they are chemically well- 


separated by many bonds. Such interactions can be evaluated by the excluded 
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volume, which is the negative integral of the Mayer f-function: 


ae a eas 


where U(r) is the potential between two monomers with a distance of r. According to 


this equation, the Mayer f-function describes the difference between the Boltzmann 
factors of two potential fields: one is U(r), the other is U(r) = 0 as r > oo. The 


excluded volume v is an integral of f(r) over the whole space, 


‘je ‘ AiGr= [ox (-<) ope [ox (-72) dr, 


Therefore, the excluded volume is a quantity that determines whether the net in- 


teraction between two monomers is attractive(v < 0) or repulsive(v > 0). 

In polymer solvent, the interaction between monomers is affected by the solvent 
molecules. When the monomers like to stay with the solvent molecules more than 
with each other, v is positive and thus this solvent is called “good” solvent. When 
the monomers are more likely to stay with each other, v is negative and thus the 
solvent is called “poor” solvent. When v = 0, the state is called “6-state” and the 
solvent is called “6-solvent”. In a polymer melt, the polymers are the 6-solvent of 
themselves, because the monomers cannot distinguish which polymer chain they 
come from. Therefore, the static properties of the polymer chains in melts can be 
evaluated by the same scaling of ideal chains. Our discussions throughout the thesis 


will be restricted to the melt state. 


1.2.4 Gaussian Chain 


Since the end-to-end vector of the ideal chain follows Gaussian distribution, one can 
also conclude that the distribution of the vector between monomers 7 and j also 
follows Gaussian, as long as the chain segment length between them is much longer 
than the Kuhn length b. It seems that the freely jointed chain model is already 


good enough to give the static properties of flexible polymers in melt state. But this 
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Figure 1.4: The discrete (a) and continuous (b) Gaussian chain model. 


model is not the simplest in mathematics. The simplest model is called Gaussian 
chain model, which assumes the chain follows Gaussian distribution over all length 
scales. 

The discrete Gaussian chain model is shown in Fig. [1.4{a). N + 1 beads are 
connected by N harmonic springs, whose potential is given by 


3kpl 


Upona(Ri+1 _ R;) = ~Op2 


(Rai —R,, 1=0,1..N 


The vector between any two monomers also follows Gaussian distribution 
3 a 3(R; — R,)? 
6-ARRD= (ara) oe (“SE ) 
This notation can be written into a continuous formula, in which R; — R;_1 is 
replaced by OR,,/On, where n ranging from 0 to 1 is the coordinate of the points on 
the continuous chain (Fig. [1.4{b)). Then the last equation can be written as 


PIR) =Cep (de [ ae (S| . (1.4) 


where C is a constant. 


1.3. Observables for Polymer Dynamics 


In this section, we will introduce a few observables to describe the relaxation dy- 


namics for linear rheology. 
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1.3.1 Stress Relaxation 


One important observable to describe the relaxation dynamics is the stress tensor 
\1’7|. To define the stress tensor of a homogeneous system, one can consider a volume 
V, which is divided by a hypothetical plane perpendicular to 6-axis. The boundary 
of the volume along (-axis is 0 and L. The stress tensor og is the component of 
force per unit area on the plane (a-axis is orthogonal to 3-axis and thus parallel to 


the hypothetical plane): 


vag = Sah = fat (Fal), 


where A is the area of the plane, and the angular bracket is the ensemble average, 


the force Fy, is force which the upper part exerts on the lower part through the 
plane. Defining X, as the integral for F,.(h), we have 


Xq= [ dh(F(h)) = eS (Fine [ dhO(h — Rmg)O(Rng — hy) 


where Finna is the component force along a-axis that monomer n exerts on monomer 
m, Rmg is the coordinate of monomer m on (-axis, h is the coordinate of the plane 
along G-axis, O() is the Heaviside step function which restricts that the hypothetical 
plane must be in the between of monomer m and n. Since the integral is non-zero 
only when the hypothetical force is in the middle of monomer m and n, the last 


equation could be written as 


Xo = © Finn (Rng — Rmg) O (Rng — Ps) 


Exchanging m and n, and using Newton’s third law: Finne = —Fnma, we have 


Xq = € > Finn (Rng — Rp) [@ (Rng — Rp) + © (Rip — rao 


1 
= (-5 Fhe (Ring — ra) 
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Defining F,,, as the sum of the forces acting on monomer m: F,, = 3°, Finn, the 


stress tensor of a homogeneous polymer melt is given by 


1 1 1 1 
Cap = AL (-5 d FrneFms + 9 > Fmt — ate (3 Fl ’ 


After a small step strain yo along x axis, in xy plane the stress relaxes in the 


form of 
Oxy(t) = YoG(t), 


where G(t) is so-called stress relaxation modulus or stress relaxation func- 
tion. In experiments, G(t) can be easily measured by rheometer with oscillating 
mode, where the oscillating strain changes with time, yt) = yo sin(wt), where w is 
the frequency of oscillation. The stress response due to elasticity is instantaneous to 
the strain, and thus the in-phase response is called storage modulus G’(w). The 
stress response due to viscosity is proportional to the shear rate, and the out-of- 
phase response is called loss modulus G’”’(w). G’(w) and G”(w) are respectively 


sine and cosine Fourier transform of G(t): 
Ct) Su | saedGbar Cy Su | cos(wt)G(t)dt. (1.5) 
0 0 


In experiments, G’(w) and G”(w) are useful to characterize the relaxation regimes. 
In computer simulation, it is easy to measure G(t). According to the fluctuation- 
dissipation theorem [72], the response of a system in thermodynamic equi- 
librium to a weak external field (e.g., magnetic and electric field) is the same as its 
response to a spontaneous fluctuation, and the change of a physical quantity is a 
linear function of the field. Using the fluctuation-dissipation theorem, G(t) of an 
isotropic system is given by 

Gt) = oF (oes lt)oas(0)), (1.6) 

B 

where the angular brackets indicates ensemble average, a and (§ are any two orthog- 


onal directions and G(t) is averaged over all pairs of a and (3. 
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1.3.2 End-to-End Vector Relaxation 


In computer simulation, the end-to-end vector relaxation function of a polymer chain 


can be easily calculated by the autocorrelation function: 
(Re(t) - Re(0)) 
(R.(0)*) 


In experiments, R, of the polymers whose molecular dipoles are parallel to backbone 


on) = 


can be directly measured by dielectric spectroscopy [8]. By applying an external 
electric field, the polymer would be polarized. For each chain, the sum of the 
molecular dipoles is equal to the induced polarization P(t), which is proportional to 


end-to-end vector R,(t). Therefore, the end-to-end vector relaxation is given by 
(P(t) - P(0)) 


2 = Bop) 


(1.7) 


1.3.3. Mean-Square Displacement 


Mean-square displacement (MSD) is a quantity to describe the diffusion of 


particles. The mean-square displacement of monomer 7 is given by 

g(i,t) = ((Ri(t) — Ri(0))’). (1.8) 
The center-of-mass mean-square displacement is given by 

ga(t) = ((Rem(#) - Rem(0))*) 


which is widely measured in experiments to give diffusion coefficient of the chain 


[73] [74]. 


1.4 Rouse Model 


The Rouse model provides analytical solutions for almost all observables in unen- 
tangled polymer melts. In this section, we will introduce the analytical solution for 
the relaxation time of the Rouse model, which is related to the discussion in chapter 


For other quantities, their scaling behaviours will be briefly introduced. 
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1.4.1 Rouse Chain 


The schematic plot of the Rouse chain is identical to the discrete Gaussian chain 
as shown in Fig. [1-4{a). Consider a chain with N + 1 beads that are connected by 
N harmonic springs of an average bond length b. The potential of the chain U, is 


the sum of the free energies of the springs 


U, (Ro... Rw) = =5- a (Rug=B: (1.9) 


In melt state, the non-bonded interaction applied to a given monomer due to the 


collision of surrounding particles could be represented by a random force satisfying 
(fi(t)) =0, (fiat) fia(t’)) = 2kpTE6ij5a96 (t — t’) . 


where 7 and 7 are the monomer index, a and £ are the Cartesian components, and € 
is the friction coefficient of the beads. Then the equation of motion of the monomer 


7 is given by 


dt? OR; 
In an “overdamped” system, the left term of the above equation can be neglected, 
giving the Langevin or stochastic equation: 


dR;  _ OU (Ro,..., Rw) 
dt OR; 


+ f,(t). (1.10) 


1.4.2 Rouse Modes 


In matrix notation, the Langevin equation for all monomers can be combined into 


a simple equation: 


N 
dR; 3kpT 
ar sear S° AugR, + f(t). (1.11) 


j=0 
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A;; is a N-order connectivity matrix, 


where only the non-zero elements are presented. It is difficult to directly give a 
analytical solution for Eq. because the motion of the monomers are correlated 
with their neighbours. But some independent behaviours could be decomposed. 
These behaviours are coherent on certain length scales. For example, the mechanical 
wave on an oscillated string is a coherent motion of the points on the string, which 
could be described by a cosine function. On a chain, similar coherent motions can 
be extracted by diagonalizing the matrix A;;. Hence the equations of motion are 


transformed to the normal coordinates 
dX 
Oa = Shp hig + ipl) PHU; (1.12) 


where X, defines the coherent motion of a chain section of length N/p and is called 


Rouse mode, which is a function of R,, 


N ; 
1 mp(t + 1/2) 
x, = —— R; —_____— ]}. 113 
P N+1 > = ( N+1 ae 
The dynamics of each mode is equivalent to the diffusion in a quadratic poten- 
tial field. The elastic coefficient of the potential field k, and the effective friction 
coefficient €, are respectively given by 


0 (p = 0) 
kp = 4 24kpT(N +1) sin? ( ug) ) hte: (1.14) 


b 
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and 
--f WDE @=0) i 
2N+1)€ (p=1...N) 
The inverse transform from X, to monomer coordinate is given by 
~ (i + 1/2) 
R; = Xo +2 5_X, cos (m= (1.16) 


p=1 

For a Rouse chain, X, is a sum of cosine functions whose wave length varies 
in respect of p (See Eq. [1.13}). Apart from Xo which represents the diffusion of 
the center of mass, other modes represent the coherent motion of a chain section 
containing N/p monomers. For example, X, describes the relaxation over the whole 


chain length. 


1.4.3 Relaxation Times 


Eq. is a stochastic differential equation (SDE) of a Ornstein-Uhlenbeck 


process, 


€,dX, = —k,Xpdt + \/2kpTE,dW, 


where T’ is the temperature, and W is a Wiener process. The solution of this 


SDE is 


X,(¢) = X,(0) exp (—t/T,) + ./2Dp i exp (-* — “| dW’, 
P 
Te sl Mess DD, = “ee (1.17) 


The autocorrelation function of X, is 


(Xpa(t)Xqa(t!)) = bpq5ag (Xpa(0)Xpa(0)) exp (-! - ") 


t—t t+t 
+DpTp (exw (-—") — exp (- = )) (1.18) 
Tp Tp 
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Let t' = 0 and t > 0, it becomes 


(Xpa(t)Xqa(0)) = Spq5ais (Xpa(0)Xya(0)) exp (-=) » CQ) 2te = 7 
(1.19) 


Thus 7, is the characteristic time or relaxation time of Rouse mode p. According to 


Kgs. {1.14} and |1.17} 7, is given by 


2 ee Tp £b? N? 
— me : 1.20 
 2kepT O(N +1)) 302 kpTp? vey) 
The relaxation time for the slowest mode is the so-called Rouse time, 
fb? 1 &b? N2 
=7,= rm) : 1.21 
SE kel NON SED) Ore kal eel) 


For the fastest mode, the relaxation time is 


oF. i aN Eb? 
= —_ ~) : 1,22 
™N [kp O(N 41) / ~ 12keT ee) 
In the qualitative discussion on scaling behaviours [75], it is also common to use 79 


TR 
N?’ 


which is in fact 2.5 times smaller than Ty. 


7 (1.23) 


1.4.4 Monomer Mean-Square Displacement 


When p = 0, the solution of Eq. gives the center of mass mean-square displace- 


ment of the chain: 


Dromelt) = ((Xolt) — Xo(0))*) = Greys 


This equation simply follows the Kinstein-Smoluchowski relation, thus the whole 
chain could be visualized as an effective Brownian particle with a diffusion coefficient 
D&kpT/NE. 

The monomer mean-square displacement, gi(t), is more complicated. At the 


time scales t < 7, the monomers are not aware of that they are part of a chain, thus 
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experiencing free diffusion. After 7), the monomers follow Rouse motion. According 


to Eqs. and |1.23} the relaxation time for a Rouse mode 7, can be written as 


ie N\? 
One) a 
Pp Pp 


During 7,, the monomers on average move a distance of the order of (N/p)b*, which 
is the mean-square size of the sections containing N/p monomers. Substituting t for 


Tp, we have 


p= (my =N a (1.24) 


t t 
which means, at the time scale t, the index of relaxed modes must be higher than 


p. The monomer mean-square displacement is of the order of the mean-square size 


of the section involved in the coherent motion at this time scale: 


git) = (=) . b? (t <t <TR). (1.25) 


TO 
For t > Tr, the monomer follows free diffusion, the mean-square displacement is 


linear in time. In summary, g;(t) for Rouse model has 3 regimes: 


bs t< TO 
gilt) _ fe, T <t<7R (1.26) 
t, t>TR 


1.4.5 Stress Relaxation and Viscosity 


After a step-strain, each unrelaxed mode contributes the energy of the order of kgT 
to the stress relaxation modulus. Therefore, at the time scales 7) < t < Tp, the 
stress relaxation modulus G(t) is proportional to the number density of sections 


with N/p monomers at time 7,, where p is given by Eq. |1.24| leading to 


db db t —1/2 


where @ is the polymer volume fraction, G(t) decays with the time as t-!/?. The 
analytical solution of G(t) for a long chain is given by 


kpT 2t 
GG)= oe exp (-=) 


p=1 
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According to this equation, we can estimate the power laws in the short and terminal 
regimes. At the time scales shorter than 79, the fastest mode has not yet fully 
relaxed. According to the expansion of single exponential, e~?/ = 1 — t/(279), 
the early relaxation shows linear decay, G(t) ~ 1 — t/(27)). Beyond the longest 
relaxation time Tr, the stress relaxation modulus has a single-exponential decay: 


kpT 2 
G(t) ¥ ee exp (-) (Ste) 


In summary, there are 3 regimes in the stress relaxation function of Rouse chain: 


t 
1- — t< 
QTn’ mn 
G(t) ~ Te ™7 <t< Tr, (1.28) 


1 
N exp(—2t/7R), t > TR 


The viscosity of the Rouse chain can be calculated from the time integral of G(t): 


a i PTR 4 1 
ica [ CO)e = a ODE. = N 


which explains the linear relationship between viscosity 7 and molecular weight MW 


for unentangled polymers. 


1.4.6 End-to-End Vector Relaxation 
The analytical solution of the End-to-end vector correlation function is derived using 
the inverse transform of Rouse mode (see Eq. (1.16), and the end-to-end vector is 
N 
Re =Ry—Ro=—4 S~ cos (sy) x 
p=1,odd 
Substituting this equation into Eq. the correlation function becomes 


b(t) = vary & tan? (orn) exp (-=} 


ons odd 


CHAPTER 1. INTRODUCTION Zi 


Since eee tan~? (mp/2(N + 1)) = N(N + 1)/2, this function would start from 


1. ®(t) also has 3 regimes: 


6kpT 

i t t< 
ENB? 2 
4 i 


From t = 0 to 70, ®(t) decays from 1 to 1—1/(2N), this early time decay is hard to 
observe when N > 1. Therefore, unlike the stress relaxation modulus G(t), ®(t) is 


not sensitive to early time relaxation. 


1.5 Tube Model 


In this part, we will introduce the tube model. This model originated from the 
study of rubber elasticity [15], then was expanded to uncrosslinked systems by 
de Gennes [16], and formulated by Doi and Edwards [17]. Afterwards, theories for 
entangled polymers have been for half century primarily based on the tube model, 


which is considered as one of the most successful models in polymer dynamics. 


1.5.1 Mean-field Tube 


The diffusion of an entangled polymer is much slower than that of the free Rouse 
chain due to the topological interactions with surrounding chains which prevent 
them from crossing each other. Consider a chain confined in a polymer network, the 
topological constraints imposed on the probe chain could be described by a set of 
obstacles scattered around the chain, as shown in Fig. {1.5} These obstacles do not 
affect the static properties of the probe chain, but significantly change dynamics. 
Since the transverse fluctuation of the chain is restricted to a “tube-like” region 
(dashed curves in Fig. [1-5), the chain can be assumed as confined in an effective 


tube with a constant diameter. Such mean-field assumption captures the essence of 
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Figure 1.5: A schematic plot of the tube model. 


entanglement and allows relatively simple mathematical treatment. Along the tube 
contour, the shortest path between the two chain ends is called the “primitive path”, 
as shown by the red curve in Fig. On the primitive path, the projected chain 
is called the “primitive chain”, which has a contour length L. 


Defining the primitive chain in a continuous manner, say R(s,t) is the position 


of the chain segment s at time t, s € [0,..., Z]. The tangent vector on s is given by 
OR(s,t 
u(s,t) = on) ) 
Os 


The original tube mode has two assumptions about the primitive chain. The first one 
is that the contour length of the primitive chain does not fluctuate. This assumption 
significantly simplifies the theoretical treatment, but introduces a non-trivial error 
(see Sec. |1.5.5). The second assumption is that the correlation of the tangent 
vectors u(s,t) and u(s’,t) decays quickly with |s — s’|, such that the primitive chain 
is Gaussian. With these two assumptions, the primitive path could be visualized as 
a random walk with a step length a, which is the length of the tube statistical 
segments (or “tube segments”). Since the mean-square end-to-end vector of the 


primitive chain is equal to that of the 3D chain, we have 
(R2) = Za? = Nb, 


where Z is the number of tube segments or entanglements per chain. The contour 
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length L of the primitive chain is then given by 
_ Nb’ 


a 


L=Za 


According to Eq. the free energy of a chain reaches its minimum when the 
end-to-end distance is zero, which seems to imply a chain confined in a tube should 
collapse into a coil. To understand the paradoxical result, one must distinguish the 
effective tube from the infinite long tube. In the tube model, the chain ends are 
not confined in the tube, but can freely explore the polymer network. Therefore, 
the chain ends have higher degree of freedom than middle monomers, acting as a 
hypothetical tensile force applied at the chain ends to keep a contour length L: 
According to this equation, the virtual force only depends on the tube segment 
length a. 

In a tube segment, the average number of monomers is N, = N/Z, thus a is 


given by 
a= V/N.0 


In the discussion of scaling behaviours, a is also regarded as the tube diameter. 
At the length scales smaller than a, the monomers are not aware of the tube con- 
finement, and thus follow Rouse motion. Therefore, the entanglement strand with 


Ne. monomers relaxes with a characteristic time Te: 
2 
To = Ns 


which is the Rouse time of a tube segment. 


1.5.2 Reptation 


The breakthrough that makes the tube model applicable to describe polymer dy- 
namics was brought by de Gennes [16], who proposed a “reptation” model to depict 
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OR 4X0 ORS 


Figure 1.6: Reptation of the primitive chain. 


the “snake-like” motion of the chain in a fixed network. In the “reptation” picture 
(see Fig. |1.6), the primitive chain moves forward and backward along the tube 
contour with a friction proportional to the number of monomers, N€, such that the 
diffusion coefficient of the primitive chain is 


kpT 
D2, 
INE 
During reptation, the chain ends will extend out of the original tube and create 
new tube segments with random orientations. In this way, the original tube will be 
forgotten gradually. The time for the chain to completely diffuse out of the original 
tube is the “reptation time”, which could be estimated by 
(DL?) €N22 (/N\°  EN?0? 
Trep © x = 
P' De kpT \N.o 


i Vig (1.30) 


As shown in Fig. [1.6{b) and (c), a tube segment is released when either chain 
ends passes through it. The survival probability of the segment s after time t, w(s, t), 


can be obtained by solving the diffusion equation: 


OYp(s,t) _ , OY(s,t) 
ot a Os?” 


with the boundary conditions: 


~(s, 0) =1; W(0, t) = 0; w(L,t) =). 


The solution of the equation is 
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The analytical solution of 7,.) is smaller than the scaling approximation (Eq. {1.30} 
by a factor of 1/7. The fraction of the tube segments that live longer than t is an 


integral of w(s,t) over the contour length: 


1 & a. 1 a 
w== | u(s,t)ds= = > Sexp(-2 (1.32) 
L 0 mw? p=1,odd p Trep 


The tube model assumes that once a tube segment is deleted, any deformation 
associated to it is forgotten, therefore p(t) is proportional to the end-to-end vector 


relaxation function ®(t): 


1.5.3 Stress Relaxation and Viscosity 


At the time scales 7 < t < 7, the monomers follows Rouse motion, such that the 


stress relaxation modulus Gt) is the same as Eq. 


G(t) = Go (<) 1 bkpl (<) _ (T) <t < 7%) 


TO b3 To 


On the time scale t = 7, each entanglement strand contributes order of kgT to the 


stress relaxation modulus G: 


where p is the density of the melt, R is the idea gas constant. At time scales t > Te, 
the chains are relaxed by reptation. Thus the stress relaxation modulus G(t) is 


proportional to the survival fraction of the tube segments (Eq. |1.32): 


G(t) = G.u(t) = =G. > So (-2 =| ; (Loe) 


p=1,odd Trep 
Due to the factor 1/p? in Eq. |1.33} the longest relaxation time 7,.) dominates G(t). 
At the timescales T. < t < Trep, G(t) shows a plateau regime. Thus, G is also 


termed as the “plateau modulus”, analogous to a permanent polymer network. 
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The viscosity predicted by the reptation model is given by the time integral of 
G(t): 
0° ve 
No = / G(t)dt = Tp Geter: 
Since G, is independent of the molecular weight M and Tyep is proportional to M?, 
the viscosity 9 is proportional to M?. However, experimental results indicate a 
power law my) ~ M*?* [4H7], whose exponent 3.4 is significantly larger than the 


prediction of reptation picture. 


1.5.4 Monomer Mean-Square Displacement 


At time scales smaller than 7,, the monomers are not aware of the tube confinement, 


thus the monomer mean-square displacement follows the Rouse behaviour same as 
Eq. 
1/2 
t 
g(t) = (=| (ti << ee) 


TO 
After 7, the monomer displacement perpendicular to the primitive path on the 
length scales larger than a is prohibited by the effective tube confinement, while the 
diffusion along the tube contour is free. At the time scales Tr, < t < Tr, the monomers 
are involved in the coherent Rouse motion of the chain sections containing (t/ 19)! . 
monomers, thus the 1D curvilinear mean-square displacement of the monomers is of 


the order of the mean-square size of the chain sections: 


((As(t))?) © (L)" via (4) Tt <TH: (1.34) 


However, the displacement of such subdiffusive motion is much shorter in 3D space, 
because the primitive path is effectively a random walk in 3D space with a step 
length equal to the tube segment length a. Since the primitive length of the random 
walk along the tube is ,/((As (jy), the monomer mean-square displacement in 3D 


is the product of the step length and the primitive length: 


g(t) © ay/ ((As (t))?) & @ (=) m mol <r 
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At t = Tr, gi(t) & Nb? indicates the monomer mean-square displacement is of the 
order of the mean-square end-to-end distance R?. At the time scales larger than 
Tr, the monomers are involved in the reptation of the primitive chain, thus the 1D 


mean-square displacement is 


((As(t))?) » Det = = Ti REX Tess: 


The corresponding mean-square displacement in 3D space is 


gi(t) © ay/((As (ay = at Ti Tee 


At the time scales t > Trep, the whole chain follows free diffusion, g(t) ~ t. 
In summary, the mean-square displacement of the monomers g;(t) have 4 regimes 


in reptation model: 
ie T9 <t< Te 
ge Te <U< TR 
gi(t) ~ me (1,35) 
es ee ee 


iy: bo Tp 


1.5.5 Contour Length Fluctuation 


In the reptation model, the primitive chain or primitive path is assumed to have a 
constant length L. Doi suggested that the contour length can fluctuate due to 
the Rouse motion of the chain. Thus, one can expect that the characteristic time of 
the contour length fluctuation (CLF) is the Rouse time 7g. According to Eq. 
the inherent curvilinear leads to contour length fluctuation: 

t 


1/4 
) (te <0 << mR) 
TO 


1/2 
(1L(#) — L(O)P)"? = 6 
The contour length fluctuations at the early times slightly reduce the stress 


relaxation modulus G(t), because a fraction of tube is released: 


(L) — noe LO) 1- ae (<) < (Get < Te): 


1/2 


G(t) = G. 2G 
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Since (tg/Te)'/* is proportional to N/?, the difference between G. and G(tR) de- 
creases as the chain length N increases. Therefore, the contour length fluctuation 
is more important on relatively short chains. 


According to Eq. {1.36| the stress relaxation modulus at the Rouse time of the 


Ni fay N, 
(ey ~G.(1—pr/—), 
7 (2) PN 


where p is a coefficient of order unity [75]. The reptation time is also reduced by a 
2 
factor of @ — pb N./N) 


2 2 
(L?) Ne Wet. Ne 

Trp @ = [1 —- —| 27 - —]|. 
rep Dz. bh N 0 N. bt N 


The viscosity can be estimated by the product of the dominate reptation time and 


chain is lower than G,: 


G(TR) x ee 


the stress relaxation modulus at Tyep: 


3 
tokpT N3 N. 
No G (tren) Teen x a W2 1- LA/ 


Over a reasonable range of molecular weight, this equation predicts that the molar 
mass dependence of viscosity, ny ~ M34 , which explains the difference between 
the predictions of reptation model and the experimental results. For very long 
chains, the contribution of the contour length fluctuations become less important, 


thus the exponent will slowly fall back to 3.0 with increasing chain length. 


1.5.6 Constraint Release 


The model combining reptation and contour length fluctuation can describe the 
dynamics of polymer chains in a matrix of much longer chains or in a fixed network. 
In a monodisperse system, the chains imposing topological constraints to a probe 
chain also experience reptation and contour length fluctuations in their own tubes. 


Therefore, the tube segments of the probe chain can fluctuate when some of these 
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Figure 1.7: The schematic plot of constraint release. 


chains move away by CLF or reptation. As shown in Fig. [1.7{a), the departure of 
a neighboring chain is represented by the disappearing of one obstacle (red), which 
leaves a vacant volume of dimension of a?. The probe chain can explore this volume 
as shown in Fig. [L.7{b). If a new chain moves in, as shown by the inserted obstacle 
(blue) in Fig. [1.7{c), the probe chain effectively hop a distance of a. The local 
rearrangement of the tube segments due to the exchange of the surrounding chains 
is called constraint release. 

The constraint release process could be modelled by the Rouse motion of the 
tube consisting of N/N, segments, called constraint release Rouse (CR Rouse) mo- 
tion [25]. The hopping frequency of tube segments is inversely proportional to the 


reptation time 75,, of the surrounding chains that impose constraints [76] [77]. There- 


ep 


fore, the relaxation time of CR Rouse motion of the tube is 


2 
P N 
TCR ~ Trep N. . 
e 


For monodisperse systems, Tor is proportional to N°, which is much longer than the 
reptation time, Trep ~ N°. Therefore, one should not expect the constraint release 
to affect the terminal relaxation time T;ep or the end-to-end relaxation ®(t) strongly, 
because they are dominated by reptation of the chain itself. 

Combining both reptation and constraint release, the diffusion coefficient of a 


probe chain could be approximated as 


2 2 
pi Ate! 4. Sel (1.37) 
Trep TCR 
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In monodisperse systems, T;rep is much shorter than 7cR, and thus the second term 
on the right-hand side of Eq. is negligible. In polydisperse systems, the con- 
tribution of constraint release on chain diffusion can be non-trivial. Consider the 
binary blend of polymers with drastically different molecular weights, where the long 
chains are confined in the matrix of the short chains. When the reptation time of 
the surrounding chains is much shorter than the reptation time of the long chain, 
the diffusion coefficient is dominated by constraint release effect. 

Due to the broad distribution of the hopping rates of the tube segments, signif- 
icant amount of constraint release events happen much earlier than the reptation 
time of the probe chain, whose effect needs to be taken into account for predict- 
ing stress relaxation, even for monodisperse systems [79]. The stress relaxation 


caused by the CR Rouse motion of the tube is proportional to t~!/? (see Eq. |1.27): 


Ger(t) ~ (= i 


Accordingly, certain fraction of the stress would be released at T¥,,,. 


The CR contri- 
bution to the stress relaxation modulus of monodisperse systems was first deduced 


by Graessley [80]: 
G(t) = Gault)R(), 


where G(t) is proportional to the product of the fraction of the surviving tube at tf, 
u(t), and contribution due to CR Rouse motion of the tube, R(t). In the “double 
reptation” model, R(t) ~ u(t), we have G(t) ~ p(t)? [20]. Rubinstein and Colby 
presented a self-consistent model to describe the stress relaxation of the binary 
blends, in which the stress modulus is calculated using a linear mixing rule of the 


stress relaxation moduli of the two components: 


G(t) = Ge (drui(t) Ri(t) + dsus(t)Rs(t)) , 


where “L” and “S” represent the long and short chains, and ¢ 1,5 refer to their volume 


fractions. 
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1.6 Multiscale Computer Simulations 


1.6.1 Molecular Dynamics Simulation 


Due to the rapid advances of the computer technology, molecular dynamics (MD) 
simulation is widely used in the study of polymer dynamics. With well-defined 
potential field, the motion of the molecules can be obtained by integrating the 
equations of motion, e.g., the Newton’s equation. Comparing to coarse-grained 
models, MD simulations based on bead-spring models carry much more fine details 
and can provide microscopic understanding that are generally not accessible by 
experiments. However, MD simulation requires much longer simulation time than 
the coarse-grained models. A typical MD simulation system for polymer dynamics 
usually has to simulate more than 104 ~ 10° particles, with a simulation time span 
over nine decades. Such computational cost strongly limits their applications to the 
entanglement dynamics at large time scales. 

For studying dynamics in polymer melts, a well-developed generic bead spring 
model is the Kremer-Grest (KG) model [49] [50]. In this model, the non-bonded 


interactions are given by the truncated-shifted Lennard-Jones (LJ) potential, 


g\ 12 a\6 1 
ioe |(F) =e) a ole (1.38) 


0 CT. 
where r is the distance between two beads, € is the depth of the potential well, o 
is the diameter of the beads, r. = 2!/8o is a typical cutoff-distance beyond which 
the potential is zero. € and o are used as the units for potential and distance. The 
corresponding LJ time unit is t,j = o(m/e)'/?, where m is the mass of the beads. 
In the KG model, the bonding potential between two adjacent particles in a chain 


is given by the finitely extensible nonlinear elastic (FENE) potential, 


1 mye 
—-kRiln (: — (x) r= Ao 
Urenn(7) = - Ro (1.39) 


oe) r > Ro 
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Entanglement 


Slip-link 
Slip-spring 


Anchor point 


Figure 1.8: The schematic plot of the slip-spring model. 


where k = 30¢/o? is the spring constant, Ry = 1.50 is the maximum bond length. 
FENE potential can effectively prevent chain crossing, thus the KG model is appli- 
cable to entangled polymers. In the KG model, the particle density is p = 0.8507, 
and the average bond length is (py? = 0.970. 

In addition, using harmonic bending potential is considered as an effective way 
to introduce chain stiffness and consequently more entanglements to the chain. A 


typical bending potential is given by 
Unend(@) = ko (1 — COS 0) 5 
where the coefficient kg is usually set to be 2€ or 3e. @ is the angle of the neighboring 


bond vectors. 


1.6.2 Slip-Spring Model 
Model Construction 


The basic building block of the slip-spring model is the Rouse chain. Consider a 


Rouse chain constituted by N + 1 monomers. The bonding potential between the 
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neighbouring monomers 7 and 7 + 1 is harmonic: 


3kpl 
b? 


Uren t,¢ - 1S (R; — Rig) ‘ 


The topological constraints imposed by the surrounding chains are represented by a 
set of virtual springs, as shown in Fig. |1.8} Each virtual spring has one end (“anchor 
point”) anchored in space, and the other end connected to the chain via the slip-link. 
The slip-links are the small rings that the chain can pass though, due to reptation 


or CLF. The potential of virtual springs is harmonic, 


3kpl 


= We (a;-—Rz,)’, 


Uss(j) 


where a; is the coordinate of the anchor point of the j-th virtual spring, x; is 
the index of the bead that the slip-link sits on, N, is the number of monomers in 
the virtual spring. The average number of slip-springs on each chain is equal to 
Z°S = N/N®*. There are thus on average one slip-link per chain segment of N5° 
monomers. N5° and N, are a pair of adjustable parameters in the slip-spring model, 
their combination determines the plateau modulus G®* and terminal properties [66]. 
The total energy of a chain is a sum of the potential of the bonds and the virtual 


springs: 
ZSs 


N-1 
Oasis = S° Usend(t,t 7 1) Tr S- Uze( 7): 
i=0 j=l 


The equation of motion of the Rouse beads in the slip-spring model is similar to Eq. 


1.10 
dR; OU chain 
ot OR; + Ce) ay 


Algorithms for Slip-Link Motion 


In the slip-spring model, an algorithm is required to govern the diffusion of the slip- 
links along the chain. There are two versions of slip-spring models with different slip- 
link motion algorithms. The original version uses Brownian dynamics that allows the 


slip-links to slide along the chain continuously [82]. The friction of slip-links 
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€, is artificially defined, but must be much smaller than the frictions of the Rouse 
beads €, i.e., &; = 0.1€ in Ref. [66]. A small friction &) ensures that the slip-links 
does not affect the chain diffusion. But in accompany with it, a small time-step At 
must be employed, e.g., the standard time-step is 0.0179 for integrating the equation 
of motion of the Rouse beads and slip-links, where 79 is the slip-spring unit time. 

In order to increase At, an updated slip-spring model employs Monte-Carlo 
method to govern the motion of the slip-links [83]. Specifically, at each time-step, a 
slip-link hops between the nearest neighbouring monomers with a frequency of f**. 
The hopping is governed by the Metropolis-Hasting algorithm: the hopping from 
monomer i to j (j = i+ 1) is accepted with a probability of exp(-AE/kpT). AE 
is the change of the potential energy due to the hopping: 


3kpl 


AE= Ne 


(R5 — R? + 2a-(R;—R,)), 


where a is the position of the anchor point. When AF < 0, the hopping will always 
be accepted because exp(—AE/kpT) > 1. 


Constraint Release 


With the assumption that the entanglements are binary events, the slip-spring model 
can easily incorporate constraint release effect by simulating an ensemble of chains. 
The algorithm is based on the “duality” of the constraint release: disentanglement 
between a pair of chains takes place when either chain moves away by reptation or 
contour length fluctuation. Particularly, each slip-link is paired with another slip- 
link sitting on another chain in the ensemble. When one slip-link is deleted from the 
end of one chain, the paired slip-link will also be released instantaneously, leading 
to a CR event. To ensure a constant number of the slip-springs in the ensemble, a 
new pair of slip-springs would be added back immediately after the removal of the 
old pair: one slip-spring is attached to the end of a randomly selected chain, while 
the paired one would be added on a random bead of another chain. Because the 


slip-links cannot sit on the same bead, the newly added slip-slinks must be added 
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on the beads without slip-links. It should be noted that the slip-links on the same 
chain are not coupled, because the fraction of the self-entanglement in this model is 


not fully consistent with the real system. 


Chapter 2 


The First-Passage Time Problem of 


One-Dimensional Rouse Chain 


2.1 Overview 


The simplest branched polymer is 3-arm symmetric stars. In the description of 
the tube model, star arms are confined in their own tubes. Reptation is highly 
suppressed because one arm has to simultaneously drag all other arms into the same 
tube to reptate. In order to relax, arms must retract some distance towards the 
branch point along the primitive path, and poke out to create new tube segments. 
By repeating this process, the old tube segments that have been visited by the arm 
free end will be forgotten, whereby the stress associated with them would be relaxed. 
An arm is fully relaxed when its free end retracts all the way back to the branch 
point. Due to the steeply growing quadratic potential field, the arm retraction 
process, essentially analogous to contour length fluctuation, is a thermally activated 
process, whose relaxation time can be obtained from the solution of the first-passage 
problem (or Kramers problem [81]) for a diffusing chain [27]. 

Consider the primitive chain as a 1D Rouse chain along the tube contour with one 


end fixed to the branch point and the other end stretched by a virtual tensile force 


42 
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\1’7|. The arm retraction problem is a multi-dimensional first-passage problem, whose 
dimensionality is equal to the number of Rouse modes. Milner and McLeish 
solved this problem by treating the whole chain as one bead attached to the branch 
point via a harmonic spring, which essentially reduces the multi-dimensional problem 
into a 1D case by only considering the slowest mode, which might overestimate the 
relaxation time. 

As the simplest stochastic model of polymer dynamics, the 1D Rouse model al- 
most has everything solved analytically, such as the stress and end-to-end relaxation 
functions. A notable exception is the first-passage problem, which in fact is a general 
challenge for the study of rare events |84]. To cope with such problems, many ad- 
vanced numerical methods have been developed to accelerate computer simulations 
in order to provide reliable numerical solutions [68H70} [85}H87|. These numerical 
methods, which have been shown to be remarkably successful in other scientific ar- 
eas, should also be applicable to the first-passage problem of the 1D Rouse chain 
model. 

At the beginning of this chapter, we introduce the solution of the Milner-McLeish 
theory without constraint release and an asymptotic solution for multi-dimensional 
first-passage problems. In order to examine these two theoretical predictions, we 
employ advanced numerical methods, including the forward flux sampling (FFS) 
and weighted ensemble (WE) method, which have been reported with excellent 
performance for solving other first-passage problems. Before being applied to the 
arm-retraction problem, both methods are tested on a simple 2D case, which is the 
simplest multi-dimensional first-passage problem. Then we choose the FFS method 
to investigate the extension problem of 1D Rouse chain, which is analogous to the 
arm-retraction problem. The mean first-passage time T(z) is measured when the 
free end of the chain extends over a certain distance z away from the origin (fixed 
end). The results show that the mean first-passage time is getting shorter if the 


Rouse chain is represented by more beads, which validates the prediction of the 
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asymptotic solution. 


2.2. Theoretical Solutions of Arm Retraction 


2.2.1 From Retraction to Extension 


Consider an arm confined in an effective tube as shown in Fig. [2.1{a). At length 
scales shorter than the tube diameter a, the monomers are not aware of the confining 
tube, and thus the dynamics follows Rouse motion. Above the length scale of a, 
one can use the primitive chain to describe the dynamics. Consider the primitive 
chain as a 1D Rouse chain along the tube contour, whose one end is fixed to the 
branch point and the other end is free to move, leading to creation or deletion of 
tube segments (see Fig. [2.1{b)). Because the arm free end is able to explore more 
directions in the polymer melt than the middle monomers, the primitive chain has 


a non-zero average length due to the entropy gain at the end: 
L=(N/N.)a 


where N is the number of monomers, a is the tube diameter, and N, is the entan- 


glement segment length. The entropic tensile force acting on the chain end is given 
as (Fig. [2.1{c)): 
3kpl— 3kpT 
fea Nb? a 


where kp is the Boltzmann constant, T' is the temperature, and b is the Kuhn length. 


Define z as the length of the tube sections that have been visited by the free end 
during arm retraction, as shown in the Fig. [2.1{b) and (c). The energy barrier for 
the arm end required to retract a distance z is given by [27]: 


_ 3keT 
a 


U(z) (2.1) 


This potential function takes the same form as for a 1D Rouse chain with one end 


fixed to extend over a distance z, as shown in Fig. [2.1{d). Therefore the Rouse 
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(a) 


(b) 


(c) 


Z 


" \ aoamaaateatal 


Figure 2.1: A schematic plot of the transformation from the arm-retraction problem 


to an extension problem of a 1D Rouse chain. 


chain extension is a reverse process of arm retraction, and so can be mathematically 
treated in the same way. Due the steeply growing energy barrier, the first-passage 


time to retract a distance z can be approximated by 


r(z) = Toexp (2) | 


where 7) was intuitively considered as the Rouse time Tr of the 1D Rouse chain with 


one end fixed [26]. 
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Figure 2.2: Rouse chain with one end fixed. 


2.2.2. Rouse Chain with One End Fixed 


Consider a Rouse chain with one end fixed, whose N beads are connected by N — 1 
harmonic springs with an elastic constant k = 3kgT/b?. The first bead 7 = 1 is 
connected to the origin by another harmonic spring, as shown in Fig. Using a 


N-order connectivity matrix, the equations of motion of beads can be written as: 
€)dR; = —kA,jR,dt + V 2kpT Ed W;,, i= 1, ‘ite (2.2) 


where & is the friction coefficient of beads, W is a Wiener process, A;; is the element 


of the connectivity matrix: 


Multiplying R; by matrix ® which consists of the eigenvectors of A, the equations 


of motion are transformed into normal modes: 


€,dX» = —kyXpdt + ./2kpTE,dW,, (2.3) 
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where the Rouse mode X, is given by 
c= mp (i + 1/2) 
x, = ———— R,; sin | ————_— ], = 1, 2, <ias IV, 2.4 
p Noi 2 sin ( N+1/2 ) : ae) 
and the friction coefficient €, and the elastic constant k, are given by 


24kpT(N 4+ 1/2) sin? Ge 1/2) ) 
2( 


(2.5) 


The inverse transform from X, to monomer coordinate is given by 


— 1/2) 
=2> XK . 2.6 
. pin (A N+1/2 ) ee) 
For each mode, the relaxation time 7, is given by 


f&  § ob". 5 { m(p—1/2) 
_ TSS) eee: 
P= ker (MON eye) pe 


Thus, the relaxation time of the fastest mode is Ty & €9b?/12kpgT, and the Rouse 
time TR is given by, 
AEgb?,.N? 


TR=7T1 ~ 372 kel’ 2.0) 


which is 4 times larger than Tg of the chain with both ends free. 


2.2.3. The Kramers Problem in Arm Retraction 


According to Eq. the extension of 1D Rouse chain can be decomposed into 
independent Rouse modes, whereby the bead coordinates R; are converted into X, 
via Eq. [2.4] (italic rather than bold font is for 1D case). Since each mode corresponds 
to one degree of freedom, this extension problem can be treated as a particle diffusing 
in an N-dimensional space. Therefore, the arm retraction problem is equivalent to 
the following first-passage problem: a particle is injected at the origin and deleted 


when it reaches the absorbing boundary satisfying 


24, sin (ae): z>a, (2.8) 
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while the potential in each dimension is given by 
1 
U(X,) = 5 heXp- (2.9) 


In the Milner-McLeish theory [37|, this multi-dimensional problem is simpli- 
fied into a 1D case by coarse-graining the whole chain into a single bead connected to 
the origin via a harmonic spring, whose elastic constant is 3kgT'/Nb?. The equation 


of motion of the bead is 


te == xt f(t), (2.10) 


where z is the coordinate of the particle, f(t) is the random force satisfying (f(t) f(t’) 
eke (t — t'), Eo = NE&q/2 is the effective friction coefficient. €.¢ is half of the 
chain friction because the center of mass travels half the distance of the free end. 
The relaxation time 7)1(z) is the first-passage time that the bead extends over a 


distance of z. 


2.2.4 Exact Solution of 1D Kramers Problem 


In this subsection, we start from a general solution of the Kramers problem without 
specifying the dimensionality. Then we will introduce the exact solution of 1D 
Kramers problem and use it to calculate Tyy1(z). 

Consider a Brownian particle in a deep potential well, U(R). To escape from 
the potential well, the particle has to overcome an extremely high energy barrier at 
the boundary R,, U(R,) >> kgT. Define W(R,t) as the probability density to find 
the particle at coordinate R and time t. w(R,t) can be obtained by solving the 


Smoluchowski equation: 


Op 


ap = V (VU(R) + keT'V) ¥. 


If there is no absorbing boundary, w is independent of ¢ and follows Boltzmann 


distribution: 


w(R) = ex (- ) N= [ox (- ) dR. (2.11) 
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With an absorbing boundary, one can find a steady state solution for w. In steady 


state, w is time-independent. The current J(R) is given by Fick’s law, 


J(R) = —7U(R)VUIR) = “2° Vu(R). (2.12) 


Assuming that ~ follows Boltzmann distribution in most places apart from the small 
region close to the absorbing boundary, it is convenient to change the unknown w(R) 


in this case to 


#(R) = o(R) exp (so) (2.13) 


Substituting ¢(R) into Eq. 2.12} we have 


In a 1D case, R is reduced to x. J is constant everywhere, which can be replaced 


J(R) = 


by J. Then we can rewrite the above equation into 


dd(a) J 
dx kpT 


U(x) 
kpT 


exp(——_). 


At the absorbing boundary x = z,, we have ¢(x,) = 0, thus the solution of the last 


equation can be written as 


The mean-first passage time is given by 


= 5 | veda = 7 foo (-73) ds 


Substituting Eq. to the last equation gives 


a oe ke U(x)\ [* Ue). 49 
T= el. dx exp (- 9 / exp (2 ) dx 


The inner integral is dominated by x close to x,, one can expand U(x) near 2,, 


U(x) + U(as) + ( — x5)U'(a5). The outer integral is dominated by the minimum 
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of the potential where « = a, we can expand U(x) © U(ay) + (U" (xp) /2) (x — a). 
Then the mean first-passage time is given by 


at QkpT a - U(x) — U(2p) 
~ T(x) \ Ua) > ( ke ) 


T 


In Milner-McLeish theory, we have U(x) = 3kgT'x?/(2Nb’), € = N&q/2, rp = 0 


and x, = z. Setting kgT to 1, we get the relaxation time: 


qe 387 
Tum (Zz) & 4/6" s exp ( 5 ) , (2.15) 


where s = z/V Nb?, TR is given by Eq. 


2.2.5 Asymptotic Solution 


Base on the early works of Kifer [88], the asymptotic solution for the N-dimensional 
Kramers problem was obtained and proved rigorously by Meerkov [84] [89]. With a 


similar idea, Cao et al. proposed an asymptotic solution, which gives 


HOVE € 2rkpT det(d) ae) . (2.16) 


~ UL(xs,¥s) Vo det(Aj;) kpT 
In this equation, the potential field U(R) is redefined by rotated coordinates, as 
shown in Fig. After rotation, the z-axis is perpendicular to the absorbing 
boundary 2, Y = {y;} are the axes of other dimensions orthogonal to x. The 
current close to 2 is assumed parallel to the x-axis, and dominated by the minimum 
of the potential on Q: (a,, Y,). A is the Hessian matrix of the potential computed at 
the minimum of the potential in the original coordinates. 4 is a covariance matrix. 
Applying it to the arm-retraction problem, the relaxation time 7(z) in the limit of 


large N is given by 


T(z) 


5/2 N 2 
eet YN? exp ( ad ) (2.17) 


7 W6N 2z 2Nb? 
where Tr is given by Eq. This asymptotic solution is 2/N times smaller than 
the prediction of Milner-McLeish theory without constraint release (Eq. |2.15). 
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Figure 2.3: Coordinate rotation according to the absorbing boundary in the asymp- 


totic theory. The left and right plots are before and after the coordinate rotation. 


2.3. Advanced Numerical Methods 


In the previous sections, we have shown some analytical solutions for the first-passage 
problems. These solutions can only predict the first-passage time in the limit of the 
infinitely high barrier. Apart from these analytical solutions, one can use numerical 
solutions, such as computer simulations, to solve the problems. However, rare events 
occur in a extremely low frequency, making them hardly observed in brute-force 
simulations. 

In the past a few decades, great efforts have been put into the development of 
numerical methods that can accelerate the computer simulations on rare events. 
Starting from the initial idea of “reactive flux” 91], several classes of methods 
have been developed, such as the “transition path sampling”, the “conformation 
dynamics”, and the “reactive trajectory sampling”. In this part, we will test two 
advanced methods that have been reported with excellent performance on first- 
passage problems, namely the forward flux sampling and weighted ensemble 
\69] [70) [85H87) methods, which belong to transition path sampling and reactive 


trajectory sampling respectively. 


CHAPTER 2. FIRST-PASSAGE PROBLEM OF 1D ROUSE CHAIN 52 


2.3.1 2D Kramers Problem 


To test their performance, we will use the FFS and WE methods to solve the simplest 
multi-dimensional Kramers problem: the escaping time of a Brownian particle from 


a 2D harmonic potential well. The potential field is given by 
din in, bp 8 
U(x, y) = sat” + By (2.18) 
The equation of motion of the particle is 
OU 
€dr = a + \/2EkpTdW, (2.19) 
r 
where r = (x,y). The absorbing boundary locates at 
22 E+ yu, (2.20) 


This toy model is intrinsically analogous to the extension model of a 1D Rouse chain 
with only 2 beads. 
The minimum of the potential on the absorbing boundary is at the point (5, ys): 
_ _ Byz _ Bez 
Ba + By’ 
When U(x, ys) >> kgT’, the current on (2, ys) will dominate the flux through absorb- 


Ls 


ing boundary. According to the asymptotic solution in Eq. |2.16} the first-passage 


time for this toy model is 


r(2) _ (% — Coe a ( bbe ) 


When the elastic coefficients are equal, 3, = 6,, this problem will be reduced into 


1D case. Thus we should set different coefficients, e.g., 8; = 1 and 8, = 10. Then 


the analytical solution is given by 


11\3? (kpTa\ 71 52? 
(2) = (=) ( - ) 5 exp (=) : (2.21) 


which will be compared with the FFS and WE results to test their accuracy. 


CHAPTER 2. FIRST-PASSAGE PROBLEM OF 1D ROUSE CHAIN 53 


2.3.2 Forward Flux Sampling Method 


Unlike most transition sampling methods, FFS requires no prior-knowledge about 
the phase space density, but needs a clear definition of the reaction coordinate. 
The reaction coordinate could be any quantity that can characterize the transition 
process, whose choice in principle only affects the efficiency. Using the reaction 
coordinate, a sequence of non-intersecting interfaces A; (i = 0...m) can be defined 
to divide the phase space into many layers. Through these interfaces, a set of 


consecutive samplings are performed instead of the brute-force simulation. 


3 


do A Nia Ay Aint 
(b1) (b2) 


Figure 2.4: (a) The schematic plot of the interface definition in FFS for a general 
transition from state A to B. (b) The schematic plot of two FFS stages. 


Consider a general transition from state A to B as shown in Fig. [2.4{a), and 
assume that the potential of state B is much higher than that of state A. The space 
between the two states have been divided by a set of interfaces, while the first 


interface Ao and the last interface ,, are the borders of A and B. The standard FFS 
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method proceeds in two stages. In the first stage, a long simulation starting from 
A is continuously performed to explore the space. During this stage, the trajectory 
will cross the interface A; many times, as shown in Fig. [2.4{b1). Among the crossing 
points, we count those last crossed interface Ag rather than \,, and recorded them in 
the database for 1, e.g., the 3 labeled crossing points in Fig [2.4(b1) will be counted. 
The result of this stage is the attempt frequency: 


Yo = No/To, (2.22) 


where No is the number of counted crossings and 7 is the simulation time. 

In the second stage, a series of short consecutive samplings are performed from 
interface A; to Am—1. For the interface \;, the sampling simulations start from 
the points on the interface (randomly selected from the database), as shown in 
Fig. [2.4{b2). These simulations finish when their trajectories either reach the next 
interface ;,; (successful run), or go back to the first interface A) (unsuccessful 
run). The first arriving points of the successful runs on A;+; are recorded in the 
database for Aj;41, e.g., points 1 and 2 in Fig. [2.4{b2). The fraction of successful 
runs gives an estimate of the probability to progress from one interface to the next, 
P(\ia|Ai) = Ni/Mi, where M; is the total number of runs from 4; to A;41, and N; 
is the number of successful runs. Thus the mean first-passage time from state A to 


a interface X,, (n >> 1) is given by 


T(An) = = I Povald) : (2.23) 


As shown in Fig. the application of FFS method to the 2D Kramers problem 
is straightforward. Defining the reaction coordinate by AX = x+y. The first and last 
interfaces are Ay = 0 and A, = z. Other interfaces \;, 2 = 1...m— 1, are equally 
inserted between Ap and A,,,, demanding that the transition from A to B must cross 


all of them in sequence. 
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Figure 2.5: Application of the FFS method to the 2D Kramers problem. 


Time-Step 


In a Brownian dynamics simulation, the continuous trajectory is represented by the 
discrete movement steps, whose variance is proportional to the time-step At. Thus 
the precision of the simulation results strongly depends on At, because the variance 
of the movement steps determines the probability of an unobservable crossing event 
as shown in Fig. [2.6{a): a movement step (solid arrow line) shows no crossing while 
its real trajectory (dashed arrow line) has actually crossed the interface. Without 
counting these crossings, the first-passage time will be surely overestimated. In 
Fig. the first-passage time 7 from direct simulation are plotted at different At 
(circles). One can find 7 decreases quickly with the reducing At. Fitting the data by 
a parabolic function, the extrapolated 7 at At = 0 is 665, which is still smaller than 
T = 690 obtained at At = 10~*. Therefore, the systematic error is non-negligible 
even with a small time-step. 

To cope with this problem, Ottinger presented an algorithm to predict the 
probability of the unobservable crossings. As shown in Fig. [2-6{a), 1, and ly are the 
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(b) 


Figure 2.6: Setting parameters according to the Ottinger’s algorithm: (1) time-step 
At, (b) interface distance AA. 


distance of the particle to the interface before and after a time-step. The crossing 


probability of the real trajectory during this step is given by 


Poros = exp(— Ns (2.24) 


where D = kgT’/€ is the short time diffusion coefficient of the reactive coordinate. 
Thus one can use a random number uniformly generated on [0..1] and compare it 
with Paoss to determine whether the crossing happened or not. The results of the 
direct simulation optimized by this algorithm are the cubics in Fig. |2.7| For all time- 
steps, the results agree well with the extrapolated 7 at t = 0 from direct simulations 


without Ottinger’s algorithm. 


Interface Distance 


Intuitively, a small interface distance is in favour to enhance the transition rates be- 
tween the neighbouring interfaces, especially for those with a steep potential barrier 
. However, the Ottinger’s algorithm requires that the interfaces distance cannot 
be too small, because it gives only the probability to cross one interface. When 
the interface are very close, the probability that the particle crosses more than one 
interfaces in At is non-negligible, which may cause a significant systematic error 


by underestimating the transition rate. To prevent this problem, a safe interface 
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Figure 2.7: A comparison between the first-passage time 7 obtained from the direct 
simulations with and without Ottinger’s algorithm. The simulations are performed 


on the 2D toy model with an absorbing boundary at z = 4 


distance should be calculated according to Poross. 
Define I, and 15 as the distance of the particle to the second interface before and 
after a time-step At, as shown in Fig. [2.6(b). The probability to cross the second 


interface P2,q must obey 


ae Ad 
Pong = exp (-23) < exp (-Fa) , (2.25) 


where A. is the interface distance. In the 2D case, the diffusion coefficient D = 2 


when the friction coefficient on each axis &) and the energy kpT’ are set to unity. 
Assuming the crossing probability of the 2nd interface is negligible when it is smaller 
than 0.01, the safe distance can be obtained at different At, e.g., AA > 0.3 when 
At = 0.01, and AX > 0.09 when At = 0.001. 

It must be noted that the distance AXg between Ag and A, should be larger than 
other interface distances. Because the crossing points on A, are obtained during 


a continuous simulation in the first stage, when Aj and A; are too close, these 
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crossing points will be strongly correlated. In order to reduce the correlation, one 
should choose a larger distance, AXg, and a longer simulation time, 7p. In our 
simulations, we set AAp to be 2.0, and other interface distances smaller: AX = 0.5 


when At = 0.01, and AX\ = 0.1 when At = 0.001. 


Systematic Error and Averaging Method 


In a multi-dimensional problem, the exploration of the phase space is very expensive 
(there could be more than 100 dimensions in the 1D Rouse chain problem!). In 
the meantime, the FFS method works in a consecutive manner, which restricts 
the application of advanced techniques, such as the parallel computing. Therefore, 
instead of getting the first-passage time by a very long FFS simulation, it is better 
to conduct a series of short independent FFS runs (on different CPUs) and calculate 
the first-passage time by averaging their results, such that the computational cost 
on each CPU is much cheaper. 

In order to investigate the errors due to the averaging methods, we perform the 
FFS simulations on the 2D case with an absorbing boundary at A,, = 14, where 
the potential barrier is higher than 88kpg7’. For all independent runs, the sampling 
number from each interfaces, M;, must be identical to ensure an equal statistical 
weight of each trajectory. For convenience, M; for all interfaces is set to be a 
constant. We compare the averaged results using two different sets of parameters: 
(1) M; = 10* and Ng, = 100, (2) M; = 10° and Ng, = 10, where Ng, is the number 
of independent FF'S runs in each case. One can find both sets have exactly the same 
total sampling numbers, Ng,M;. M; in first set is 10 times smaller than in the second 
one, thus the variance of its results is bigger, which means larger statistical error. 
For both parameter sets, the time step At is 0.001, and the interfaces distance AX 
is 0.1. 

In Fig. T(A) is normalized by a factor Aexp™!(Umin(A)), such that 7(A) 
should approach to a plateau when A > 1 (see Eq. (2.21). The triangle symbols in 
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A Arithmetic mean: M=10000 WN,,=100 
vy Arithmetic mean: M=100000 N,,=10 
o Harmonic mean: M=10000  N,,=100 
o Harmonic mean: M=100000 N,,=10 
--- Theoritical prediction 


1(A) A EXP (Umin(A)) 


Figure 2.8: First-passage time 7(A) obtained by arithmetic and harmonic means 


from the 2D Kramers problem. 


Fig. [2.8]represent the results by arithmetic mean, 


where 7;(A) is the mean first-passage time obtained in k-th independent run. When 
A < 4, 7(A) is much higher than the asymptotic solution, due to the low energy 
barrier at the early stages. With increasing , T(A) is suppose to gradually converge 
to the predicted value (dashed line in Fig. [2.8). However, after \ = 8, the normalized 
T(A) values obtained by using the first parameter set deviate from the plateau value 
and increases sharply. For the second parameter set, similar deviation happens at 
slightly larger A (after 4 = 10), and the growth rate is lower than that using the 
first set. Also, such unexpected behaviour seems related to the sampling number 
M,. If we keep increasing M;, the deviation of 7(\) from predicted value happens 
at larger . 


We consider this deviation as a systematic error arising from the averaging 


CHAPTER 2. FIRST-PASSAGE PROBLEM OF 1D ROUSE CHAIN 60 


method. Arithmetic mean is valid only when these independent runs have exactly 
the same statistical weights. In FFS, there are a few samples which have very long 
T(A) (many times of the expected value), which will dominate the arithmetic mean 
and so lead to unrealistic large mean value. We can compare this averaging method 
with one direct simulation. In a continuous long run, 7(A) can be calculated by the 
total simulation time divided by the number of first crossings on the boundary 
(those last crossed Ap now crosses A), because the time spent from back to Xo is 
negligible comparing to 7(A). In a continuous run, the trajectories that arrive \ with 
a shorter period of time will also return to the Ay immediately to start another trip to 
A, which means these faster trajectories have higher probability or higher statistical 


weights. It implies that 7(A) should be calculated by the harmonic mean 


= Nas 
le). 


As shown in Fig. |2.21} the circles and squares representing the harmonic means are 


T(A) (2.26) 


very close to each other and both converge to the asymptotic line. Additionally, a 
general observable A associated to the interface (eg. the first-passage point on the 


interface) could be averaged by the following equation: 


ON Ag) rad) 
AQ) = SNe tra) ae 


where A;(A) is the observable obtained in k-th FFS run. Replacing A,z(A) by 7% (A), 
this equation is exactly the same as Eq. 

It should be noted that the harmonic mean is not an universal averaging method 
for all first-passage time problems. For example, we will use arithmetic mean in 
the first-passage problem of slip-spring model discussed in next chapter. This is 
because there are many starting states, while the transitions between them are also 
rare events. In that case, the faster arriving trajectories have the same statistical 


weight as the slower ones. 


CHAPTER 2. FIRST-PASSAGE PROBLEM OF 1D ROUSE CHAIN 61 
2.3.3 Weighted Ensemble Method 


The WE method was first proposed by Hurber and Kim [69], who employed it in 
Brownian dynamics simulations to study binding process in protein. Rojnuckarin 
et al. [85] used this method to explore the configuration space of folded coarse- 
grained proteins. Later, Zhang et al [70] extend the WE method to a broad class 
of stochastic dynamics. Recently, Darve et al. improved the resampling 
algorithm, making it applicable to various transition problems. 

In the WE method, the phase space is divided into many non-overlapping do- 
mains, as shown in Fig. [2.9] These domains could be hexagons (Fig. [2.9{a)), layers 
(Fig. [2.9{b)), or any other arbitrary shapes. For convenience, the domains in the 2D 
Kramers problem are defined by the non-intersecting interfaces similar to the FFS 
method; therefore, we call the domains as “layers” in discussion. In simple words, the 
WE method works as running a lot of parallel brute-force simulations with different 


priorities. Such priority is weighted by so-called “resampling” algorithm. 
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Figure 2.9: Domains in WE method. 


Before applying the WE method to a specific problem, it is necessary to introduce 
the general idea of resampling. Resampling is a flexible algorithm that duplicates 
or kills the sample elements without changing their distributions. Consider a col- 
lection of numbers which follow the Gaussian distribution with an average of 0. For 


convenience, we set the weights of the numbers to be 1. If we randomly delete half 
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of the numbers, and double the weights of the rest (or keep their weights to be 1 but 
duplicates the rest of the numbers), the distribution and the sum of weights does 
not change. Supposing we are more interested in the numbers larger than 0, another 
resampling algorithm with two steps could be taken: (1) find all the numbers smaller 
than 0, delete half of them and increase the weights of the other half to 2; (2) find 
all the numbers larger than 0, randomly duplicate half of them and let their copies 
take half of their weights. With these resampling steps, a new collection of numbers 
are created, which is dense in the particular part of interest and sparse on the rests. 
In the WE method, the mean first-passage time to an absorbing boundary 4 is 

calculated by 
(A) == (2.28) 


where ¢ is the simulation time, N) is the total number of trajectories which have 
crossed the interface X, w? is the weight of j-th trajectory that has crossed A, the 
constant Wo is the total weight of the trajectories in the system. Similar to Eq. 


2.27, a general observable associated to A, A(A), can be calculated by 


Ay = (2.29) 


Resampling Algorithm 


Apparently, one can design his resampling procedure. However, it is rather diffi- 
cult to find a balance between efficiency and accuracy. In the early WE method 
for Brownian dynamics, new trajectories are created by splitting the weights of the 
older ones, leading to a strong polarisation on the weights of trajectories in each 
layer, i.e. some trajectories have very large weights while the others are very small. 
Such polarisation makes the results hard to converge to the correct value. Recently, 
an advanced resampling algorithm without suffering such polarisation has been pro- 
posed [87]. Briefly, this resampling algorithm tend to delete the elements with small 
weights and split the elements with large weights; therefore, this method not only 
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keeps a constant number of elements but also keeps the weight of the elements equal 
in each domain. 

We use Fig. [2.10{a1) and [2.10{a2) to explain the resampling algorithm adopted 
in the 2D case. In Fig. [2.10{a1), there are 3 trajectories labelled by different colours 
in the first layer A,. Define the weight of the red, blue and grey trajectories are, 
respectively, wr, Wp and we (wr > Wp > Wg). Suppose the expected number of 
trajectories in each layer is My, = 2, the first step of the resampling algorithm is to 
kill one trajectory and give its weight to another. Because the red trajectory carries 
the greatest weight, the choice would be made between the other two trajectories. 
Since the survival probability of a trajectory is proportional to its weight, one can 
uniformly generate a random number on (0...1]; if it is larger than w,/(wp + We), we 
kill the grey trajectory and add its weight to the blue one, and vice versa. In the 
example of the layer A, in Fig. [2.10{a2), we keep the blue trajectory, whose weight 
becomes Wy =wptw,. If uy, = Wy the resampling is finished, otherwise further 
steps must be taken. For example, in layer A» of Fig. [2.10{a1), the red trajectory 
has larger weight than the blue one, w, > wp, and their average is W = w, + Wp. 
Because w,; > @, we split the red trajectory into two, their weights are w, = w 
and w, = w, — w, respectively. Now, it turns into the situation in layer A,, and 
one should repeat the previous step. The iteration terminates when the trajectory 


number is My, in layer A and their weights are equal. 


Parameters 


For the 2D case, the simplest layer definition is equally dividing the phase space by 
a set of interfaces as in the FFS method. But this definition can be further improved 
by adjusting the width of the layers. Since the total weight of trajectories in each 
layer roughly follows Boltzmann distribution at steady state, it is reasonable to set 
density of the interfaces sparse at the bottom of the potential field and dense close 


to the absorbing boundary, as shown in Fig. b). Because only the absorbing 
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Figure 2.10: (a) Resampling algorithm. (b) Definition of interfaces for the WE 
method in the 2D Kramers problem. 


boundary A,, requires the judgement of crossing by the Otinger’s algorithm, small 
interface distance is not an issue in WE simulations. In our case, the coordinates of 


interfaces are given by 
eae FH 1-9: (2.30) 
m 


which is a parabolic function whose maximum locates at A», = z. 

Apart from the layer definition, the performance of the WE methods also depends 
on a few parameters, i.e., the resampling frequency f,,, the number of layers m, 
and the expected trajectory number ma, in each layer. We separately investigated 
these parameters on the 2D Kramers problem with an absorbing boundary at A = 
9. As shown in Fig. 2.11{a), higher resampling frequency is always in favour to 
raise the hitting rate on the absorbing boundary, N)/tepu, where tepy is the CPU 
time. Therefore, the resampling procedure should be taken every time-step. In Fig. 
[2.11{b), we fix M, to 20 for all layers, and change m from 10 to 30. The real-time 
first-passage time T“° has been normalized by the asymptotic value r**” from Eq. 


such that it should converge to 1. On can find that the convergence rate with 
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Figure 2.11: WE method performance with different parameters: (a) the resampling 
frequency f,s;, (b) the number of layers m, (c) the expected trajectory number in 


each layer My. 


more layers is higher than that with less layers. In Fig. 2.11{c), the number of 
layers is fixed, M, changes from 10 to 30. It is found that the convergence rate for 
My, = 30 is much higher than the other two. In addition, such convergent behaviour 
implies that one should let the system to relax before collecting the data, which can 


significantly reduce the noise at the beginning. 


2.3.4 A Comparison Between FFS and WE Methods 


In this section, we will compare the results of the FFS and WE methods on the 
2D Kramers problem. The simulations are perform on the same CPU core (Intel 
Xeon E5-2620). The parameters are set as follows. The time-step At is 0.01 for 
both methods. In the FFS simulations, the interfaces distance AX is 0.5. For each 
interface, the number of sampling M; is 10°. The number of independent FFS runs 
is Ng; = 100, which is sufficient to provides good statistics. In the WE simulations, 


layer definition is given by Eq. with m equals 20. Each layer contains M, = 20 
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Figure 2.12: First-passage time obtained from the FFS and WE methods (symbols) 
and the theoretical predictions of Eq. (dashed line). 


trajectories. All WE simulations must have a relaxation period over 600 seconds 
before collecting the data. 

The mean first-passage time 7(\) obtained by different methods are shown in Fig. 
Owing to their mechanisms, the FFS method manages to get the spectrum of 
T(A) for all interfaces in a single run, while the WE method has to set up independent 
runs for each interface A. Generally, the results from both methods are consistent 
apart from the early region 4 < 4 and the late region A > 10. When X < 4, 
the energy barrier is relatively low. Since the WE simulation is close to brute- 
force simulation, the results from the WE method are more precise for lower energy 
barrier, while the FFS method is precise only when the energy barrier is much higher 
than kg7’. When A > 10, the WE data show a sharp deviation from the plateau 
and fails to converge even after a long run. The deviations are random and could be 
alleviated by increasing m and M,. Nevertheless, it reveals a risk of the WE method 
on solving the first-passage time problems, because the users cannot predict if the 


parameters are still safe for current potential barrier. In contrast, the precision of 
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Figure 2.13: Time-cost of the FFS simulation to reach each interface in the 2D 


Kramers problem. 


FFS method is controllable by monitoring the conditional probability P(\j+1|);). 
When P(Aj;+1|;) drops to a dangerous level, one can accordingly increase M;. As a 
consequence, the FFS method is more advanced on the accuracy and stability. 

On efficiency, both methods have good performance. The WE method approxi- 
mately takes less than 600 seconds to converge as shown in Fig. while a single 
FFS simulation is even faster, which requires merely about 50 seconds to reach 
A = 12, as shown in Fig. On parallel computing, WE is stronger than FFS, 
but the latter allows a remedy by simulating lots of independent runs. Consider- 
ing its advantage on exploring the whole spectrum in one run, FF'S is still a better 
choice. 

Apart from 7(A), there are some other observables related to the interfaces could 
be obtained via Eq. [2.27]and Eq. For example, the coordinates of the average 
first arriving points on interfaces, (xfp(A), Yep(A)). Since vpp(A) + Yep(A) = A, we only 
plot yep(A) in Fig. It is found that ys)(A) obtained by FFS and WE method 
are consistent. The first arriving points are all on one side of y,(A), but gradually 


converge to it. 
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Figure 2.14: First arriving point on each interface y)(A) obtained from the WE 
and FFS simulations for the 2D Kramers problem (circles) and the minimum of the 


potential (dashed line). 


In conclusion, the FFS method performs better than the WE method in this 2D 
case, and will be employed in the following study on the 1D Rouse chain model. Nev- 
ertheless, it must be noted that WE method still has some remarkable advantages 
over FFS method, such as the flexibility on defining the domains and its mecha- 
nism close to the brute-force simulation, which makes it widely applicable to some 


complicated rare events, such as the protein folding. 


2.4 Computer Simulation Study on 1D Rouse Chain 
Model 


In Sec. and |2.2.2| we have introduced the relationship between the arm re- 
traction and the extension of the 1D Rouse chain with one end fixed. By coarse- 
graining, Milner and McLeish [80} [37] reduced the multi-dimensional first-passage 


problem into a 1D Kramers problem, and proposed an approximate solution (see 
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Sec. and |2.2.4). Cao et al. [84] proposed an exact asymptotic solution for 


multi-dimensional first-passage problems and predicted a result smaller than the 
prediction of Milner-McLeish theory (see Sec. [2.2.5). In this section, both the di- 
rect (brute-force) and FFS simulations will be applied to the 1D Rouse chain model 
to examine the analytical solutions. 

We perform computer simulations with variable number of beads representing 
the Rouse chain. The bead friction ¢9, energy kpT7’, and statistical segment length b 
are set to be unity in the simulations without loss of generality, whereby the units 
for length, time and energy are respectively b, Cob?/kgT', and kgT. The predictor 
corrector method [94] was employed in simulations. The detection of trajectories 


crossings on interfaces is improved by Ottinger’s algorithm [92]. 


2.4.1 Direct Simulation 


Direct simulation results are plotted in Fig. The horizontal axis is s = z/ VNB, 
where z is the end-to-end distance or extension length, N increases from 1 to 128. 
In a continuous simulation, when the free end last crossing so = 0 reaches s > 0 for 
the first time, its time cost is recorded as the first-passage time for s. Fig. [2.15{a) 
shows decimal logarithm of the mean first-passage time. Clearly, the time grows 
very fast with s, approximately as exp (3s?/2) as expected. The direct simulation 
can approach T(s) 10° or so. Fig. [2.15{b) shows the same data but normalized by 
T(s)8T,' exp (—3s?/2). For clarity, both predictions by Milner-McLeish theory and 
asymptotic solution are divided by the trivial factor exp (—3s?/2). Such normaliza- 
tion brings all data with one decade in vertical scale and allows clear comparison 
between theories and simulations. In particular, the direct simulation are signifi- 
cantly faster than the Milner-McLeish prediction (dotted-dashed line) when N is 


larger. 
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Figure 2.15: (a) The decimal logarithm of first-passage time 7(s) for FFS simulations 
(dots) and direct simulations (solid lines). (b) Normalized 7(s) versus s for FFS 
simulations (circles) and direct simulations (solid lines), the dashed lines are the 


prediction of Eq. Milner-McLeish Theory is shown by the red dotted-dashed 


line 
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Figure 2.16: Applicaiton of FFS method onto 1D Rouse chain extension model 


2.4.2 EFFS Simulation 


In order to extend simulation results to longer times and facilitate detailed theory 
verification and calibration, we also performed FFS simulation of the same model. 
First of all, we need to define the reaction coordinates and the non-intersecting 
interfaces, which is quite straightforward in this model. As shown in Fig. the 
extension length z is employed as the reaction coordinate. The original interface Xo 


is define at z = 0, other interfaces are placed according to 
= (1025 « G=1))N 7), 4=1,2,.:. (2.31) 


Such interface definition avoids the systematic errors due to very small interface 
distance and large statistical errors due to large interface distance. 

The simulation then proceeds in two stages. In the first stage, we run one long 
simulation for time 7p and count the number of crossings, No, of the first interface 1 
by the trajectories which last crossed interface Xo, rather than A,. Besides counting 
the crossings, we store the full chain configurations at the moments of these crossings. 

In the second stage, we run many short consecutive simulations for interfaces 1 to 
n — 1 in sequence. For the interface \;, the simulation starts from the stored points 
on the interface A; (selected at random from the database) and finish when they 


either reach the next interface \;+1(successful run), or go back to the 0-th interface 
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(unsuccessful run). The fraction of successful runs N;/M; gives an estimate of the 
probability to progress from one layer to the next, P(Aj+i1|\i), where M; is total 
number of runs from layer 7, and N; is the number of successful runs. Thus, the 
mean first-passage time is given by Eq. 

The value of MM; has a decisive effect on the statistical error of the final outcome, 
with the best strategy to increase M; for higher energy barriers between the layers 
to ensure an approximately constant N;. A simple way to determine M; is to run a 
few simulations with smaller M; and get the rough ratio of P(\;4,|A;), and estimate 
M,; for an expected N;. Ref. recommends selecting interface distances such that 
P(\i41|A;) > 0.3. Our selection satisfies this criteria. By running a quick simulation 
for N = 1, the proper M; can be obtained. Using the same /; and the same distance 
defined by s for N = 1, a proper ratio P(Aj41|A;) for larger N is also guaranteed 
since the P(A;41|A;) increases with larger N. 

In Sec. the difference between the harmonic and arithmetic means has 
been discussed. In Fig. we compare the two averaging methods on the chain 
with 32 beads. The averaging are performed on two samples: (1) N; = 10° and 
Nes = 100, (2) N; = 10* and Ng; = 10. Despite the variance ((r(s) — (r(s)))*) for 
N, = 10* is much larger than that for N; = 10°, their harmonic means are consistent, 
allowing us to reduce the computational cost of a single run by performing a lot of 
independent runs on different CPUs. This method is very useful when the extension 
length is long. 

The mean first-passage time 7(s) for different chain lengths N are presented 
in Fig. Comparing with direct simulations, a disagreement can be found at 
s < 1.5. In this region, the first-passage time given by the FFS method is inaccurate 
since the energy barrier is lower than 3.5kg7’. In the region of s > 1.5, two simulation 
methods are consistent with each other. FFS method is able to predict the first- 
passage time till s = 5.5, with the chain length up to N = 128. 

In the normalized plot in Fig. [2.15{b), all curves show a fast decay for the 
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Figure 2.17: A comparison between arithmetic and harmonic mean for averaging 


independent FFS runs of 1D Rouse chain model. 


intermediate values of s and then gradually saturate around certain transition length 
s_ with clear plateau reached in the systems with small N values. The slopes of the 
curves from the peak to s; increases with increasing N. In the mean time, the 
transition length s; also increases. One finds that the result differs from the Milner- 
McLeish theory significantly. When increasing N, the first-passage time becomes 
much shorter than their prediction, leading to the difference of a factor of 10 at s = 3 
and N = 128, and even bigger for larger s and N. This shows conclusively that 
the one mode assumption of the Milner-McLeish theory is inadequate and better 
theory must be developed. Note that this discrepancy is much bigger than the 20% 
reported by Vega et.al. [95]. 

The results of FFS simulations verifies asymptotic solution in the limit of large 
extension (z ~ Nb), where the prefactor of T(z) has z~' scaling. In the intermediate 
regime, FFS results show different scaling behaviour, which is roughly tT(z) ~ 27°. 
In Ref. [34], Cao et al. combined the asymptotic solution at large s and the one in 


the intermediate regime (calculated by so-called “minimum action path’), and 
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proposed an empirical expression: 


r(s) = (a™ | a TR EXP (=) (2.32) 


go 
where C) (NV) is given by 


327 5 . 9 1 
C\(N) = gl sin (aac): (2.33) 


C2(N) is a fitting parameter. By fixing C,(V), Co(V) is obtained by fitting the FFS 
simulation results, whereby a combination of the theory and FFS simulation leads 


to a simple expression for the first-passage time 


a a rs ae 


2.5 Conclusions 


In this work, we have studied the first-passage problem of the 1D Rouse chains, 
as a proxy for dynamics of arm retraction of isolated star polymers in a network. 
In the widely known Milner-McLeish theory, star arms are represented by Rouse 
chains inside their confining tubes and further replaced by one bead attached to 
the branch point by a harmonic spring [80]. The mean disengagement time of a 
tube segment is T(z) ~ z~'exp(U(z)/kgT). In order to check the validity of the 
Milner-McLeish theory, Cao et al. [34] proposed an asymptotic solution to solve 
the multi-dimensional Kramers problem. This asymptotic solution is only valid in 
the limit of very large extensions z ~ Nb, corresponding to a fully extended chain. 
The results show that the mean first-passage time drops significantly if the arm is 
represented by Rouse chain with more beads instead of a single bead. 

Because the large deviations of the 1D Rouse chain rarely happen, the verification 
of the asymptotic theory by a direct simulation is practically impossible. To cope 
with it, we can employ the advanced numerical methods for first-passage problems. 


Two advanced methods, i.e., forward flux sampling and weighted ensemble methods, 
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were tested on a 2D Kramers problem, which is a simplest multi-dimensional first- 
passage problem. Considering their performance on all aspects, such as accuracy, 
efficency and stability, the FFS method was chosen to investigate the 1D Rouse 
chain model. The results of the FFS simulations are in good agreement with the 
asymptotic solutions at very large extension. In the intermediate regime, T(z) shows 


richer scaling behaviours. 


Chapter 3 


Arm Retraction Dynamics of 
Entangled Star Polymers: The 
First-Passage Problem in Slip-Spring 
Model 


3.1 overview 


Development of quantitative theories for predicting the dynamic and rheological 
properties of entangled branched polymers is of both fundamental and practical 
importance. In the past decades, theoretical efforts have been primarily based on 
the concept of tube model originally proposed by de Gennes, Doi and Edwards 
[24]. Different from entangled linear polymers where reptation, contour 
length fluctuations (CLF) and constraint release (CR) are the main relaxation mech- 
anisms, reptation in branched polymers is strongly suppressed due to the effectively 
localized branch points. In the simplest case of symmetric star polymers, the stress 
relaxation is conjectured to proceed via CLF or arm retraction by which the free 


end of an arm retracts inward along the primitive path to escape from the original 
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tube segments and pokes out again to explore new tube. Since arm retraction is 
entropically unfavorable and so thermally activated, this process can be formulated 
into a first-passage (FP) problem or Kramers problem. 

A star arm retracting in a fixed network experiences a potential barrier theoreti- 
cally described by a quadratic function U(s) = vkgT'Zs? where kg is the Boltzmann 
constant, Z = M/M, is the number of entanglements per arm, M the arm molec- 
ular weight, M, the entanglement molecular weight and v a constant [26]. The 
fractional coordinate s measures the retraction depth of the arm free end. Pear- 
son and Helfand predicted an exponential dependence of the arm terminal relax- 
ation time, Tg, and correspondingly the viscosity, 79, on the arm molecular weight, 
No ~ Ta ~ exp(vM/M,) [27]. This prediction however shows large discrepancy from 
experimental data obtained in star polymer melts due to the neglect of CR effects. 
Ball and McLeish took into account the CR effects by applying the dynamic 
tube dilution (DTD) hypothesis where the relaxed arm segments are considered 
to work as effective solvent for the unrelaxed materials. Milner and McLeish further 
improved this theory by including the contributions of fast Rouse fluctuations at 
early times and solving the first-passage problem of a diffusing end monomer to 
retract a fractional distance s to get the arm relaxation spectrum 7T(s) at late times 
[30] [37]. The Milner-McLeish theory predicts the stress relaxation of symmetric 
star polymer melts reasonably well, but not the dielectric or arm end-to-end vec- 
tor relaxation function. It also encounters difficulty in using a single set of model 
parameters to describe the rheological behaviors of asymmetric star polymers with 
different short arm lengths [36]. In recent years computational models based on the 
framework of Milner-McLeish theory have been developed for describing the linear 
viscoelasticity of branched polymers with arbitrary architectures and their general 
mixtures [39] [42] [44] [47| [96]. These models have been shown to provide predic- 
tions in reasonably good agreement with experimental data for a variety of systems, 


but are facing problems in describing the linear rheology of some simple mixtures, 
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such as the star-linear blends, especially at low fractions of star polymers [44] [97]. 
Therefore more quantitative theories that can simultaneously predict different dy- 
namic and rheological properties of entangled branched polymers are still highly 
desired. The development of such theories requires the analytical solution of the 
multi-dimensional FP problem of arm retraction [34]. 

On the other hand, the coarse-grained slip-link or slip-spring (SS) simulation 
models have demonstrated strong potential in describing dynamics and rheology 
of entangled polymers [59] (64) |66] 98H10I]. For example, the single-chain slip- 
spring model developed by Likhtman can provide simulation results on multiple 
experimentally measurable observables, such as neutron spin echo, linear rheology, 
dielectric relaxation and diffusion. Using a limited number of fitting parameters, 
the predictions of this model match the results obtained from both experiments and 
molecular dynamics (MD) simulations on linear and symmetric star polymers very 
well [101]. The SS model serves as an intermediate between tube theory 
and MD simulations. As a discrete model, it not only naturally builds in all the 
relaxation mechanisms of the tube model, but also carries more system details, 
such as explicit polymer chains and entanglements [102]. At higher level of coarse- 
graining, the SS model is significantly more efficient than MD simulations using 
bead-spring polymer model, which is of great advantage in the study of branched 
polymers. Furthermore, the slip-spring model can separate the contributions from 
different relaxation mechanisms by enabling some of them while disabling others. 
This is particularly helpful for examining assumptions made in current theoretical 
models and providing valuable information for developing more quantitative models. 
One typical application is to evaluate the magnitude of constraint release effects by 
comparing simulation results obtained from entangled polymer systems with and 
without CR. 

Since deep arm retractions are rare events due to the high energy barrier, the 


time and length scales accessible to standard slip-spring simulations are still much 
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shorter than those in well-entangled experimental systems where the tube models 
are supposed to work best. Similar problems have also been seen in brute force 
simulations of many other rare events, such as crystal nucleation [104], biological 
switches and protein folding [105]. The required computational time may take 
up to several decades [106]. Advanced numerical techniques, such as the umbrella 
sampling and transition path sampling methods, have to be employed to 
accelerate the simulations. Recently the forwards flux sampling (FFS) method has 
been proposed [68} [93} [109] and proven to be successful in molecular dynamics and 
Monte Carlo (MC) studies of rare events [110]. 

In this chapter, we will combine the FFS method with the slip-spring model for 
studying the dynamics of entangled symmetric star polymers. This is a proof-of- 
concept work. To our knowledge there seems no other reported work in the literature 
on applying the transition path sampling methods to study entangled polymer dy- 
namics, especially on arm retraction dynamics. We will mainly focus on the systems 
without constraint release for the following reasons: 1) It is relatively convenient 
to implement the FFS method and find an appropriate reaction coordinate in the 
non-CR systems; 2) The terminal relaxation times in the systems without CR are 
much longer than those with CR, allowing us to test the computational efficiency 
and limit of the combined method; 3) Reliable simulation data on the FP times of 
arm retractions without CR are highly desired for examining analytical solutions 
of the multiple-dimensional Kramers problem [34]; 4) The extension of the method 
developed in the non-CR case to the CR case is fairly straightforward, as will be 
shown in Sec. With an optimized selection of the reaction coordinate, which 
is the index of the monomer that the innermost slip-link sits on, we first validate 
the proposed simulation method by producing simulation results on the terminal 
relaxation times Tg of mildly entangled star arms up to 8 entanglements in good 
agreement with those obtained from SS model simulations. The FFS simulations 


are then extended to longer arms with lengths up to 16 entanglements and so reach 
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Ta values about 6 decades beyond that accessible by brute force simulations (from 
6 x 10° to 3 x 107 SS unit time). The FP times of other original slip-links along 
the arm can be calculated using similar FF'S simulations as for the innermost one, 
which consequently provides the entire arm relaxation spectrum 7(s). Moreover, 
we propose a numerical route to construct the arm end-to-end vector correlation 
functions, @(t), and stress relaxation functions, G(t), from the discrete data stored 
at each interface during the FFS runs. Such time correlation functions are generally 
not discussed in other FFS studies. Our simulation results will contribute to the de- 
velopment of theoretical models for describing the dynamics of entangled branched 
polymers and also the general first-passage problems in multi-dimensional systems. 
The simulation methodology developed in this work should also be applicable to the 
study of rare events in other scientific areas. 

The rest of this chapter is organized as follows. In Sec. we introduce the 
single-chain slip-spring model for entangled star polymers in the absence of CR. The 
detailed description of the combined FFS and SS model is given in Sec. The 
simulation results obtained in the non-CR systems are presented and discussed in 
Sec. including the terminal relaxation times Tg, the arm retraction spectra T(s) 
and the numerical route for constructing ®(t) and G(t). In Sec. the simulation 
method is extended to study the arm retraction dynamics of star polymers in the 


presence of CR. We draw conclusions in Sec. 


3.2 Slip-spring Model for Entangled Symmetric Star 
Polymers 


3.2.1 Model Description 


In the single-chain slip-spring model for entangled symmetric stars, each star arm is 


represented by a Rouse chain with N + 1 monomers linked by N harmonic springs 
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\14| \66], as shown in Fig. One end monomer with index 0 of the chain is 
treated as the branch point which is fixed in space, while the other end with index 
N moves freely. The topological constraints on the arm are modelled by a set of 
virtual springs each of NS° beads. Each virtual spring has one end connected to the 
Rouse chain by a slip-link that can slide along the chain, and the other end, called 
anchor point, is fixed in space. The slip-spring model effectively assumes a binary 
picture of entanglements, which is qualitatively supported by recent MD simulation 
studies [52] [54] [11]. There is on average one slip-spring every NS° monomers. The 
values of NS° and N°* are adjustable for describing the intensity of entanglements. 
It should be noted that N®S* is not necessarily equal to the entanglement length N,. 
used in tube theory. Their relation will be discussed in Sec. To be consistent 
with previous publications 34} [66] [82], we choose NS° = 4 and N° = 0.5. Other 
parameters, such as the bead friction coefficient ¢), the average bond length b of the 
Rouse chain, the temperature kg7 and consequently the time scale T = €ob?/kpT, 
are all set to unity. 

The Hamiltonian of the SS model is determined by the potential energies of both 
the harmonic bonds of the Rouse chains and the virtual springs. The trajectories of 
the Rouse monomers are obtained by solving their Langevin equations of motion. In 
the original slip-spring model [66} |8T] [82], the slip-links are assumed to travel contin- 
uously along the straight lines between adjacent monomers and so can sit anywhere 
on the chain. In a latter version of this model [83], the slip-links move discretely 
by hopping from one monomer to one of its nearest neighbors with the acceptance 
rate controlled by a Metropolis Monte Carlo algorithm. The long-time behavior of 
the system is not sensitive to the details of the slip-link motion. For simplicity and 
computational efficiency, we employ the discrete motion approach in the current 
work. The slip-links are not allowed to sit on or pass through the branch points 
of the star arms. In the systems without constraint release, such as star polymers 


in a fixed polymer network, the destruction and creation of slip-links can only take 
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place at the free ends of the star arms. Different from the systems with CR |66], the 
slip-links are not coupled with each other. In addition, the slip-links on the same 
arm are not allowed to pass over each other or occupy the same monomer. This as- 
sumption introduces an effective excluded volume interaction between the slip-links, 
which is consistent with the low swapping rate between neighboring entanglements 
as revealed in a recent MD simulation of symmetric star polymer melts : 

The previous slip-spring simulations were typically carried out in an ensemble of 
chains and the total number of slip-links in the system is kept constant [66]. In the 
non-CR case, when one slip-link is deleted from a chain end, another slip-link will 
be added to the end of a randomly selected chain in the ensemble. For convenient 
installation of the FFS method, we modify the SS model for the non-CR case by 
simulating each entangled arm individually. The destruction of slip-links on a given 
arm is still incurred by the retraction of the arm free end (monomer index N), but 
the addition of new slip-links to the same arm end is now determined by a probability 


P.aq which satisfies the detailed balance condition 


(1 = pa) (Pada + PaPn-an) = Pa (Prose + 1 = pa)Pww-i) ; (3.1) 


where ps = 1/N°* is the average number of slip-links sitting on each monomer. P, ; 
is the transition probability for a slip-link to move from monomer 7 to monomer 7 and 
Poss 18 the probability for a slip-spring sitting on the arm free end to be destructed 
after one time step, respectively. Eq. thus represents the balance between the 
flux of slip-links to and from the end monomer. Assuming Py—i,v = Py,n—1 without 
loss of generality, Eq. gives Piaa © 0.167 for the system parameters NSS = 4 and 
Poss = 0.5. The modified SS model is validated by studying the static properties of 


the simulation system. 
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Figure 3.1: Slip-spring model simulation results (solid circles) and predictions of Eq. 
3.2| (open squares) on the probability distribution of number of slip-links per arm, 


P(Na, N), for symmetric star polymers with arm length N = 24. 


3.2.2 Static Properties 


The static property of the slip-spring model system of entangled symmetric star 
polymers can be well characterized by the distribution of slip-links along the star 
arms. Considering the effective excluded volume interactions between the slip-links, 
the problem is similar to one-dimensional real gas in equilibrium. The probability 


distribution of finding N, slip-links on a star arm of N monomers is simply given 


by 
PNG NS Copa) (3.2) 
here CN! ad Fig. |3.1| shows th d t bet th 
where = A 12. Es SNOWS e OO agreemen etweelh e 
N ~ NA(N-Ngl © pave ae 


prediction of Eq. and the SS model simulation results on P(N, N) for the 
system with N = 24. It can be seen that the peak value of N, is located at 
Ng = 6 in consistence with the expected average number of slip-links per arm, 
(Na) = psiN = 6. 


When there are N,) slip-links on a given arm, the probability to find the i-th 
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slip-link on the monomer « is 
i-1 WNai-i 
Gn 1G ns 


P(a,i, Na, N) _ CNet 
N 


where the numerator is a product of the possibilities to find 2 — 1 slip-links on the 
arm segment from monomer 1 to x — 1 and to find Ny — 7 slip-links on another 
segment from monomer z+ 1 to N. It should be noted that in the star polymer 
systems without CR the slip-links do not change their ordering along the star arms. 
In Eq. [8.3] the index i is considered to increase from 1 for the innermost slip-link to 
higher values toward the arm free end. Combining Eqs. [3.2] and we obtain the 
ensemble-averaged probability to find the i-th slip-link on the monomer 2: 


Paa,N)S >) Pat Nay NP Oa). (3.4) 


Na=1 

Fig. [3.2] presents the SS simulation results on P(x,i,N) for the slip-links with 
indices 7 = 1 to 6 on star arms of length N = 24, together with the predictions 
of Eq. The good agreement between the two sets of data indicates that the 
simulation systems are in equilibrium state and the randomly assigned locations of 
the anchor points can well preserve the equilibrium distribution of the slip-links. 
This is also reflected by the fact that the average number of slip-links found on each 


individual monomer is equal to p., = 0.25, see the horizontal line in Fig. [3.2 


3.3. Combined FFS and SS Method For Entangled 
Star Polymers without CR 


In the systems without CR, the topological constraints or entanglements imposed 
on a target arm are released hierarchically by the retraction of the arm free end. 
The terminal relaxation time Tg of the system is defined as the average first-passage 
time that takes the free end of an arm to reach the branch point starting from a 


random initial conformation. For well-entangled star arms, Tq grows exponentially 
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Figure 3.2: Slip-spring model simulation results (symbols) and predictions of Eq. 
[3.4] (lines) on the probabilities of finding 7-th slip-link on monomer x, P(z,i, N), 
for the symmetric star polymers with arm length N = 24. The horizontal dashed 
line shows the simulation results on the average number of slip-links found on each 


individual monomer. 


with the number of entanglements per arm Z [27]. However, full arm retraction 
rarely happens at large Z and so is generally not accessible by standard brute force 
simulations. There is also no exact analytical solution of this multi-dimensional FP 
problem. Therefore the forward flux sampling method introduced in Ref. [68] is 
employed in order to study these rare events. A successful application of the FFS 
method on studying the FP time of 1D Rouse chain with one fixed end can be found 
in Ref. [84]. 


3.3.1 Forward Flux Sampling Method 


In FFS the phase space is divided by a sequence of no-crossing interfaces denoted 
by A; (« = 0,...,7m), as sketched in Fig. [3.3(a). The starting states of the dynamic 


process are on the first interface Ao, and the reactive or terminal states are on the 
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(a) 


Figure 3.3: (a) Schematic diagram of the FFS method. The continuous yellow 
trajectory is the continuous simulation in the first stage, and the blue trajectories are 
the successful shooting simulations in the second stage; (b) Algorithm for building 
continuous arm relaxation pathways from the piecewise shooting trajectories shown 


in (a). 


last interface A,,,. These interfaces are defined by a reaction coordinate, which can 
be any parameter evolving during the process, but different choices could result 
in significantly different performance. More detailed discussion about the reaction 
coordinate is given in Sec. [3.3.2] 

The FFS method is operated in two stages. In the first stage, a very long 
continuous simulation is performed in order to calculate the frequency uo at which 
the trajectory crosses the interfaces Aj and A, in sequence. In the second stage, a 
set of consecutive shooting simulations are carried out from interface A; to interface 
Aizi for i = 1,...,m—1, which provide the transition probabilities P(\j41|A;) that 
a system starting from ; will first reach \;.; rather than return to Ao. The first- 


passage time 7, for the system starting from the first interface Ag and ending on the 
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interface A, (1 <n <m), is then given by 
_ 1 
bo Ty Pi /Ai) 


3.3.2 Reaction Coordinate 


A key issue in applying the FFS method is the choice of the reaction coordinate. 
Starting from a random initial configuration, the relaxation of a star arm in the 
system without CR proceeds by the retraction of the arm free end along the primitive 
path, passing through all the original slip-links on the arm sequentially until none left 
between it and the branch point. The terminal relaxation time is determined by the 
moment at which the innermost slip-link is released. During this process, the number 
of surviving original slip-links, Nj, on the arm drops with time from its initial value 
to 0, making it an intuitively simple choice for the reaction coordinate. Considering 
that the value of Ng is statistically proportional to the length of the surviving tube 
or primitive path, this choice would be consistent with a recent FFS study on the FP 
time for the free end of a 1D Rouse chain to reach a certain distance z from the fixed 
end where z was selected as the reactive coordinate [34]. The 1D Rouse chain study is 
closely related to the current work, because arm extension is essentially the reverse 
process of arm retraction. However, when using N, as the reaction coordinate, 
our FFS simulation results on the terminal arm retraction times are found to be 
significantly smaller than those obtained from standard SS model simulations. The 
problem arises from the difficulty in choosing equivalent starting states for the FFS 
runs. In the slip-spring model system, both the instantaneous number of slip-links 
and their distribution along the arm are subject to strong fluctuations, especially 
on the outer arm segments which undergo fast Rouse motion. In the FFS runs 
using N,) as the reaction coordinate, the starting states are collected in the first- 
stage continuous simulation as the configurations where the number of slip-links on 
the arm is equal to the ensemble-averaged value of (Na) = Nps. Shooting from 


these starting configurations, only the samples in which the values of N.; decrease 
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Figure 3.4: Application of FFS method for studying the retraction dynamics of an 
entangled star arm described by the slip-spring model. The cross (Monomer 0) on 
the left represents the branch point that is fixed in space. The interfaces \; (vertical 


lines) used in the FFS simulations are placed on the monomers of the arm. 


monotonically are considered to reach interface A; successfully. This biased strategy 
is thus in favor of the samples where the initial slip-link densities on the outer arm 
segments are higher than p,,, because in such cases the probability to lose slip-links 
at short times is higher than to gain ones. Therefore a relatively large proportion 
of slip-links on a sample arm are released by shallow arm retractions at early times, 
leaving fewer than average number of slip-links on the surviving segments of the 
primitive path. As a consequence the terminal relaxation times obtained from the 
FFS simulations are shorter than those obtained from standard SS simulations where 
the ensemble-averaged initial distribution of slip-links is uniform. These results 
imply that the reaction coordinate should be selected close to the branch point in 
order to minimize the influence of the fast fluctuating arm end. 

Since the terminal arm relaxation time is determined by the release of the in- 
nermost slip-link from the arm free end, one can track the motion of this particular 
slip-link along the arm by defining the index of the monomer that it sits on as the 


reaction coordinate. As shown in Fig. [3.4] where the 3D Rouse chain is sketched as 
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a straight line for convenience of discussion, the first interface Aj) used in FFS is set 
on monomer a (2 in this case) where the innermost slip-link originally sits on. Any 
initial configuration of the confined arm in which the innermost slip-link locates on 
monomer a@ can be taken as the starting state of the FFS simulation. The second 
last interface A,,_; is placed on the outermost monomer N of the arm, and the last 
interface A,,, is right outside of the arm free end, marking the final or reactive state 
that the arm free end has passed through the innermost slip-link and the arm is 
fully relaxed. The other m — 2 interfaces are placed on the monomers in between a 
and N. 

According to the standard FFS method, a database containing a large number of 
configurations is accumulated on each interface. In the first stage of the continuous 
simulation, the database on , is a collection of configurations whose innermost slip- 
link lastly crossed Ao before crossing A;. In the second stage, consecutive shooting 
simulations are performed from interface A; to Aj41, 27 = 1,...,m— 1 using starting 
configurations randomly selected from the database on A;. Among the M; shooting 
samples, the ones whose innermost slip-links reach \;,; before going back to Ao are 


considered as successful samples and will be stored in the database of Aj41. 


3.3.3 Simulation Details 


Apart from the reaction coordinate, the performance of the FFS algorithm can also 
be affected by some other factors. One factor is that the configurations saved in 
the database of interface \, during the first-stage continuous simulation could be 
strongly correlated with each other due to the limited running time at this stage in 
comparison with Tg. This may introduce systematic errors in the simulation results 
if the size of the database is fixed. This problem can be resolved by increasing the 
interval |; between the interfaces Ag and A;, as shown in Fig. and recording 
configurations on A, at a lower frequency w. For example, rather than recording 


every event that the innermost slip-link crosses A; when coming from Ag, one can 
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record once for every 1/w crossings. Another factor is the choices of the interface 
interval J. between A; and Aj41 (¢ = 1,...,m—2) and the number of shooting samples 
M; from each A; which determine the performance of the FFS in the second stage. 
Since J, controls the transition probabilities P(\;.1|;), a smaller l, is normally 
preferred for accelerating the shooting simulations. The number M; can then be 
chosen according to P(A;+1|A;) and the desired accuracy. 

In the current work we take 1; = 2 and lz = 1 which separate the first two 
interfaces Ay and A; by one bead and then set one interface on every bead along the 
arm. The recording frequency w has to be reduced for longer arms in order to reduce 
the conformational correlations on A, and is empirically taken to be w = 1/(N — 15) 
for arm length N > 16. Since the reaction coordinate is defined by the location 
of the innermost original slip-link, the transition probability P(A;1|\;) increases 
with 2 towards the arm free end. In order to achieve good statistics for the first few 
interfaces close to the branch point, M; should be large enough. A number of samples 
M;, = 40,000 is thus used for \;, 7 = 1,2,...,m— 1 in all of the FFS simulation 
runs. As shown in Fig. there is a non-negligible fraction of initial configurations 
where the innermost slip-links are many monomers away from the branch point 
and could be released by shallow arm retractions. The terminal relaxation times 
of such arms are thus much shorter than those of the arms with uniform slip-link 
distributions. Actually their terminal times have been reached in the first-stage 
continuous simulations without going into the second stage of FFS. These 7g data 
are still counted for calculating the distribution and the mean value of the terminal 


relaxation times. 
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Figure 3.5: Simulation results on the terminal arm retraction time Tg obtained from 


FFS and direct shooting simulations as a function of arm length N. 


3.4 Results and Discussions for Systems Without 


Constraint Release 


3.4.1 Terminal Time of Arm Retraction 


The terminal time Tq of the arm retraction process is the main and most straight- 
forward output of the FFS simulations. Fig. [8.5] presents the FFS results on Tq as 
a function of the arm length N. For comparison, we have also included the tq data 
obtained from the so-called direct shooting simulations which start from the first in- 
terface Ao and stop at the last interface \,,, without intermediate steps. These runs 
are equivalent to the slip-spring simulations using initial configurations randomly 
picked from the database on interface Ay and running continuously until the inner- 
most original slip-spring being deleted by the arm free end. For each arm length, 
the direct shooting simulation results are averaged over 10, 000 independent samples, 
while in the FFS simulations Tg is averaged over 2,000 independent runs. Since in 
each FFS run there are 40,000 samples recorded on ;, the average is actually taken 


over a much bigger ensemble than that of the direct shooting runs. Considering the 
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Figure 3.6: Average computational times required for completing a single FFS and 


a single direct shooting run on a single Intel Xeon processor. 


high computational cost, the direct shooting simulations are only performed for arm 
lengths from N = 20 to 36, corresponding to about 4 to 8 entanglements per arm 
estimated with N. ~ 4.47 as discussed in Sec. In this range of N, the FFS 
and direct shoot simulation results in Fig. show very good agreement with the 
relative differences less than 5%. The combined FF and SS method and the choice 
of the reaction coordinate are thus well justified. 

Fig. [3.6]compares the average computational times required to complete a single 
direct shooting and a single FFS run on a single CPU (Intel Xeon E5-2620). The 
direct shooting simulation is faster at short arm lengths, but its computational time 
grows exponentially with N and overtakes that of the FFS when N > 32. The FFS 
method allows us to study much longer arms. For entangled star polymers with arm 
length N = 72 in the absence of CR, the terminal relaxation time is found to be 
Tq & 2.85 x 10! which is about 8 orders of magnitude longer than that of stars with 


N = 20 and is basically inaccessible to any type of direct simulations. 
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3.4.2. Comparison with Theoretical Model Predictions 


The 7q data in Fig. show a clear exponential dependence on the arm length 
N, which is expected from the Pearson-Helfand theory for star arms retracting in a 
fixed network [27]. These results can be further compared with the predictions of 
more detailed theoretical models [37]. The Milner-McLeish theory based on 
the solution of 1D Kramers problem predicts the terminal arm retraction time in 
the absence of CR as 
5/2 2 

74 N) = Tru N) = (=) : (3.6) 
where z = \/N/N, and the arm Rouse time T(N) = 4¢o.N?b?/37?kpT. The entan- 
glement molecular weight N, can be estimated by substituting the corresponding 
FFS result on 7q(V) into Eq. As shown in Fig. the obtained N, values are 
roughly independent of N, giving N. ¥ 4.98. 

Recently Cao et al. pointed out that the first-passage problem of Rouse chain 
should be treated as a multi-dimensional Kramers problem [34]. FFS simulations of 
1D Rouse chains showed that the z~! scaling in the prefactor of Tg as predicted in Eq. 
[3.6] is only valid for very large chain extensions. In the intermediate chain extension 
regime corresponding to realistic arm retraction process, a new theory based on 
the Freidlin-Wentzell theory was proposed [112], which predicts a z~° scaling in the 
prefactor of the terminal time [Eq. 60 in Ref. [34]] 


ro) = SOLA aay (22), on 


where C(NV) is a fitting parameter. For arm lengths N > 20 we can take the plateau 
value of C(.N) = 1.2 as found in the FFS simulations of 1D Rouse chains [34]. The 
N, values calculated by substituting the FF'S data on T4(V) into Eq. [3.7Jare shown in 
Fig. which increase with the increasing arm-length and approach an asymptotic 
value of N, = 4.47 that is smaller than the N, value estimated by using Eq. 
The two theoretical models thus predict qualitatively different dependence of N, on 
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Figure 3.7: Entanglement molecular weight N, calculated by substituting the FFS 
simulation results on 7, (Fig. into the theoretical predictions of Eqs. 
(squares) and [3.7| (circles) for various arm lengths. 


N, at least in the systems without CR. Since the entanglement molecular weight 
is one of the most important model input parameters for predicting the dynamics 
and rheology of entangled polymers, this N-dependent behavior apparently needs 
further investigation for developing quantitative theories. The FFS results on Tg over 
a broad range of arm lengths should work as a benchmark for examining theoretical 
models that are typically developed for well entangled polymers. 

We note that the N, values given in Fig. [3.7|are different from that obtained by 
mapping the slip-spring model simulation results on the linear viscoelastic properties 
of linear polymer melts to the tube model predictions (N. =~ 5.7) [82]. The 
difference could be related to the use of different theoretical models for the data 


fitting and also the presence of constraint release effects in the polymer melts. 
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3.4.3 Arm Relaxation Spectrum 


Apart from terminal relaxation time, the FFS method can also be applied to obtain 
the entire relaxation spectrum of the arm. This is done in a similar way as calculating 
Ta. The only difference is to set the index of the monomer that the i-th original slip- 
link sits on, instead of that of the innermost slip-link, as the reaction coordinate. 
Accordingly the first interface Ag in the FFS method is defined on the monomer 
where the 7-th slip-link originally occupied. The FP time of the 7-th slip-link is 
recorded as 7T(X) with the fractional index X = i/(N,). The simulation results on 
T(X) are plotted in Fig. [3.8] for the arm lengths 20 < N < 44. For the systems with 
N < 36, the direct shooting simulation results are also presented for comparison. 
The agreement between the FFS and direct shooting data gets improved as the 
arm free end retracts deeper along the primitive path, i. e., with the decrease of 
the slip-link index 7 and so X. This is understandable because the release of the 
outer slip-links or entanglements is dominated by the fast Rouse-like fluctuations. 
The corresponding energy barrier is relatively low such that the FFS method does 
not work well at large X. For this reason, the most reliable relaxation spectrum, 
especially for the long arms, should be constructed by combining the FP times of 
the inner slip-links as calculated by the FFS method with the FP times of the outer 
ones obtained from direct shooting simulations. One such example is shown in Fig. 
for the systems with N = 44. The complete relaxation spectrum 7(X) can be 


directly applied to test theoretical models on arm retraction dynamics. 


3.4.4 Constructing Relaxation Correlation Functions 


In experiments the dynamics and rheology of entangled polymers are generally char- 
acterized by the dielectric relaxation or chain end-to-end vector correlation function, 
®(t), and the stress relaxation function, G(t). The calculation of these observables 
usually requires the continuous trajectories of the polymers, which are however not 


naturally available in FFS simulations, because only instantaneous configurations 
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Figure 3.8: Relaxation spectrum calculated using the first-passage times of all slip- 
links for star arms with various lengths obtained by both FFS (solid symbols) and 
direct shooting (open symbols) simulations. The dashed curves are for guiding the 
eye. The parameter X = i/ (Nj) is the fractional index of the i-th slip-link along 
the arm, which increases from X = 1/ (N,)) for the innermost slip-link to 1 for the 


outermost one. 


at the hitting points on the interfaces are recorded. Here we introduce a numerical 
route to effectively link these discrete pieces of information to construct the dielec- 
tric and stress relaxation functions. The systems of entangled star polymers without 
CR are used as examples to demonstrate the application of this algorithm. 

Fig. [3.3{b) sketches the method used to build continuous arm relaxation path- 
ways from the piecewise FF'S shooting trajectories shown in Fig. [3.3[a). Considering 
two hitting points on the terminal interface X,,,, marked as A,, and B,,, there must 
be two continuous trajectories or pathways that one can track back from them to the 
first interface Ag. As shown in Fig. [3.3{b), the pathway to state A,, is constructed 
by linking the successful shooting trajectory from the hitting point A,,_,; to A,, with 
that from A,,—2 to Am—1, and so on until reaching the point A; on the interface 4. 


The linking from A, to a start point Ag is obtained from the trajectory generated 
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in the continuous simulation in the first stage of the FFS simulations. Similarly 
the pathway to the hitting point B,, can be traced back to B; on A; and then to a 
starting point By. We note that these rebuilt trajectories are different from the true 
continuous trajectories generated in standard slip-spring model simulations, but the 
ensemble-averaged pathways obtained in these two cases should be very close, as 
reflected in the consistent ®(t) results in Fig. |3.10| From computational point view, 
the rebuilding method requires the storage of all the successful shooting trajectories 
between neighboring interfaces and also large memory for data processing. This may 
limit its application to large systems such as the fine-grained bead-spring models 
widely used in molecular dynamics simulations. 

When calculating the arm relaxation correlation functions from the rebuilt tra- 
jectories, two assumptions have been made. First, when one slip-link is destroyed 
by the retracting arm free end, the primitive path segment in between its nearest 
neighboring slip-link and itself will be forgotten immediately. This assumption is 
valid for most of the slip-links due to the discrete feature of entanglements in the 
SS model. The only exception is with the tube segment between the branch point 
and the innermost slip-link where this assumption may affect the calculation of the 
relaxation functions in the terminal regime, as discussed below. The second assump- 
tion is that the FP times on each interface follow a single exponential distribution. 
This assumption has also used in solving the 1D Kramers problem and in the Doi- 
Edwards tube model without CR [17]. Since the slip-spring model is essentially a 
multidimensional problem, we perform an extra set of simulations to examine the 
validity of this assumption. A total number of 10,000 direct shooting simulations, 
all starting from exactly the same initial configuration, are carried out to mimic a 
FFS run. The FP times for the innermost slip-link to reach different monomers, 
or different interfaces in the FFS definition, are recorded. Fig. [3.9] presents the 
probability distributions, P;(t), of the FP times on three different interfaces for the 


arms of length N = 20. It can be seen that P;(t) on interfaces with higher indexes 
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Figure 3.9: Probability distributions of the first-passage times for the innermost 
slip-link to reach different monomers or different interfaces in the FFS definition 
A; along the arm as calculated by direct shooting slip-spring simulations of star 
arms of length N = 20. All of the 10,000 simulations start from the same initial 
configuration where the innermost slip-link sits on monomer 1 next to the branch 
point. The solid lines represent single exponential fit to the simulation data in each 


case. 


can be well described by the exponential function 
| t 
P(t) = — —— 3.8 
() = Sew (—+) 38) 
where 7; is the mean FP time on the interface \;. The second assumption becomes 
valid as the arm free end retracts deeply along the primitive path. 
Following Eq. [3.8] the probability that the innermost slip-link has never crossed 
the interface A; after time t is 
0 t 
P9(t) =exp(—=), (3.9) 
and the probability that it has crossed A; at least once is 


P(t) =1—exp (-=) (3.10) 


Ti 
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Therefore the probability that the trajectory starting from Ag has crossed interface 


A; but never crossed interface A;+1 is 


t 
BO = Px (t) = Pre) = — exp (-+) + exp (- 


a 


: ) (3.11) 


Ti+1 
The time correlation function of a dynamic observable, V, whose instantaneous 


values are calculated on different interfaces can be evaluated by 


(V(t)V(0)) = (By (t)Wo+ sy P(t) W; + P& (") (3.12) 


i=1 
where W; is defined as 
h; 
1 eee 
= — nye al 
W, ae Yi (3.13) 


and h; is the number of hitting points on the interface A; out of the M;_, shootings 
from A;_1, Ag is the observable value at the k-th hitting point on A; and ve is its 
value at the corresponding starting point of the trajectory on the first interface Ao. 
For the system sketched in Fig. [3.3(b), there are only 2 hitting points on the final 
interface A,, such that h,, = 2 in Eq. 


Substituting Eqs. and into Eq. |3.12| we get 


(V(t)V(0)) = (S AW;i+1 €XP (- i ) + Hn). (3.14) 


a 
i=0 oe 


where AW;i41 = W; — Wi41. The correlation function in Eq. [3.14] is expressed as a 
weighted summation of a set of exponential functions, which is consistent with the 
tube model predictions for the end-to-end vector and stress relaxation functions of 
entangled polymers in the absence of constraint release [17]. The only difference lies 
in the last term W,,, on the right hand side of Eq. [3.14] which, if being nonzero, may 
result in an unphysical plateau after the terminal relaxation time Ty. 

The problem associated with W,,, does not exist in the tube model where the tube 
is assumed to be continuous. The arm free end can thus retract continuously along 


the primitive path all the way to the branch point and so release all the memories 
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in the original tube. As a result, W,, decays to zero for all dynamic observables. 
However, in the slip-spring model, the entanglements are represented discretely by 
the slip-links. The terminal time Tq is taken to be the time when the arm free end 
passes the innermost slip-link. It implies that the memories, such as stress and 
end-to-end vector orientation, stored in the original tube segment between this slip- 
link and the branch point are not fully forgotten right after 7g, giving to a nonzero 
ensemble average value of W,,. Actually this is an intrinsic problem for all models 
using discrete description of entanglements. We will address this problem in more 
details in a future work. For the current work, we will neglect the last term in Eq. 
[3.14] by setting W,,, = 0, which is reasonable at least for the systems with very long 
arms where the contribution from the last tube segment is relatively small. This 
approximation is analogous to the so-called disentanglement relaxation mechanism 
used in the tube-based computational models where a polymer branch is considered 
to be fully relaxed by disentanglement when there is only one or few entanglements 
left on the branch [89] 44]. 

The dielectric and stress relaxation functions calculated using Eq. with 
Wi, = 0 from the rebuilt trajectories are plotted in Fig. for arm lengths up 
to N =72. For comparison, the ®(t) and G(t) results obtained from standard slip- 
spring model simulations are also included for the systems with N < 36. The two 
sets of ®(t) curves show very good agreement in the terminal regime, indicating the 
capability of Eq. in constructing the arm relaxation functions using discrete 
FFS shooting trajectories. The discrepancy at short time scales can be attributed 
to the numerical problem that the constructed ®(t) curves start from initial values 
smaller than 1. Since the W,, term in Eq. is neglected when calculating 
@(t) = (R.(t)- R.(0)) where R, is the arm end-to-end vector, the resulted value 
of (R2(0)) at the initial time is smaller than the true mean squared end-to-end 
distance of the arms which is obtained from standard slip-spring model simulations 


and used for normalizing ®(t). The so-obtained initial values of the ®(t) curves 
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Figure 3.10: (a) Arm end-to-end vector correlation function ®(t) and (b) stress 
relaxation function G(t) obtained from standard slip-spring simulations (symbols) 


and calculated using Eq. [3.14] from the rebuilt trajectories (lines), respectively. 


are ®(t = 0) = 0.80 for arm length N = 20 and 0.92 for N = 72, respectively. 
As expected, the contribution of W,,, becomes less significant with increasing arm 
length. 

The G(t) results presented in Fig. [3.10{b) are the single-arm stress autocor- 
relation functions without considering the cross correlation contributions from the 
virtual springs [113] [114]. This choice does not affect any discussions or conclusions 


in the current work, especially when there is no constraint release effect. Different 
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from the ®(t) results, the G(t) curves calculated using the rebuilt trajectories decay 
faster than those from the slip-spring simulations, which implies the existence of 
systematic errors originated from different resources. First of all, the assumption 
that the FP times follow an exponential distribution does not apply for the first few 
interfaces due to the relatively low energy barrier, as shown in Fig. |3.9| Therefore 
the constructed G(t) curves show significant difference from the standard slip-spring 
simulation results at early times. Secondly, the calculation of the stress relaxation 
function requires a very large ensemble average for achieving good statistics. In 
the molecular dynamics and slip-spring model simulations, the instantaneous stress 
tensor usually needs to be calculated at every single time step [115]. But in FFS 
simulations the number of data points on each interface is rather limited. Thirdly, 
Eq. calculates the relaxation correlation functions using the information, such 
as the arm conformations and the locations of the slip-links, carried by the hitting 
points on each interface. These hitting points are only saved from the successful 
shooting trajectories which are in favor of the arm retraction process and corre- 
spondingly the redistribution of the Rouse beads in between the slip-links. The 
biased change of the local conformations of polymer segments leads to faster relax- 
ation of the stress, because G'(t) depends on the local bond or segment reorientation. 
The end-to-end vector correlation function is less sensitive to this problem, because 


®(t) is determined by the relaxation of the vectors linking neighboring slip-links. 


3.5 Extension of the combined FFS and SS method 
to Systems with Constraint Release 


The combined FF'S and SS method can be extended to entangled polymer systems 
with CR by adjusting the definition of the reaction coordinate. In the standard slip- 
spring model [54] [66], constraint release is included by coupling the slip-links sitting 


on different polymer chains or arms into pairs to represent the binary entanglements. 
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When one slip-link is deleted from the free end of an arm, its coupled partner is also 
deleted regardless of its location, which results in a CR event. This means that 
for FFS simulations the originally innermost slip-link could not be used to define a 
reaction coordinate alone for exploring the entire arm relaxation spectrum, because 
this slip-link may be destructed by a CR event before reaching the arm free end. To 
resolve this problem, we refer to a recent slip-spring simulation work on entangled 
symmetric star polymers with CR [54]. There it was shown that the relaxation of the 
original tube segments, and correspondingly the relaxation of the arm end-to-end 
vector, is dominated by the first-passage times of the so-called tube-representative 
(TR) slip-links, which are the original slip-links finally released from the arm free 
end. The other original slip-links which are destructed from the middle of the arm 
by CR events only contribute to stress relaxation. For determining the terminal 
relaxation time of the arm end-to-end vector, we only need to find the moment 
when the last tube segment held in between the branch point and the innermost TR 
slip-link is released by the arm free end. Since it is not known in advance whether an 
original slip-link will be deleted by the arm end or by CR, we can define the reaction 
coordinate as the index of the monomer which the innermost surviving original slip- 
link sits on. In other words, if at time t the innermost original slip-link was deleted 
by CR, the reaction coordinate will be immediately shifted from the monomer it 
sat on to the monomer occupied by the nearest original slip-link, because the latter 
becomes the innermost surviving original slip-link. This procedure will continue 
until the last surviving original slip-link is destructed by the arm free end and so 
the terminal relaxation time 7,4 is reached. 

The ensemble-averaged terminal relaxation times, 74, obtained in the FFS sim- 
ulations with the modified definition of the reaction coordinate are presented in 
Fig. together with the terminal relaxation times of the arm end-to-end vector 
relaxation functions as obtained from standard slip-spring model simulations and 


the mean FP times of the innermost surviving original slip-links as obtained from 
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Figure 3.11: Simulation results on the terminal arm relaxation times Tg obtained 
from the FF (open squares) and direct shooting (open circles) simulations, together 
with the terminal times of the arm end-to-end vector correction functions calculated 
from standard slip-spring simulations (open triangles), in the systems with constraint 
release. For reference, the FFS results on 7g for the systems without CR (solid 


squares, same as in Fig. 3.5) are also plotted. 


the direct shooting simulations. The three sets of data show very good agreement 
within error bars, which effectively validates the proposed FFS method. The com- 
bined FFS and SS method can thus provide quantitative predictions on the terminal 
relaxation times of entangled star polymers either with or without CR over a broad 
range of arm lengths that are surely needed for the development of quantitative 
theories for entangled branched polymers. The construction of the relaxation cor- 
relation functions, ®(t) and G(t), in the CR cases is rather complicated and will be 


left for further study. 
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3.6 Conclusions 


We present an application of the forward flux sampling method in combination with 
the slip-spring model on studying the arm retraction dynamics of entangled star 
polymers. The single-chain slip-spring model originally developed for describing 
entangled linear polymers has been extended to model symmetric star polymers. As 
a proof of concept, we start with the systems without constraint release where the 
entanglements or slip-links can only be created on or deleted from the arm free ends, 
making the FFS method conveniently applicable. Two possible reaction coordinates 
for the FFS simulations have been tested. The choice of the index of the monomer 
that the originally innermost slip-link sits on is found to provide FFS simulation 
results on terminal relaxation times Tg in good agreement with those obtained in 
direct shooting simulations for mildly entangled stars with arm lengths up to 8 
entanglements. The FFS simulations are then performed to study the terminal 
relaxation of much longer arms (up to 16 entanglements) that are not accessible 
by any direct simulations, especially considering the exponential growth of Tq with 
the arm length in the absence of CR. The FFS results on 7g over such a broad 
range of arm lengths allow direct comparison with the predictions of theoretical 
models which are typically developed for well entangled polymers. The entanglement 
molecular weight N, extracted from such comparison is found to have an arm-length 
dependence. 

In addition to the terminal arm relaxation time, the first-passage times of all 
other original slip-links on a given arm can also be conveniently calculated by defin- 
ing the reaction coordinate as the index of the monomer that the interested slip-link 
sits on, which in turn provides the entire relaxation spectrum of the arm. For mildly 
entangled arms the FFS results on the FP times show good agreement with direct 
shooting simulation data for the deep entanglements or inner slip-links, but some 
discrepancy exists for the shallow ones, because the FFS method does not work 


well at low energy barriers. The reliable relaxation spectrum of long star arms 
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thus should be constructed by combining the FP times of the inner slip-links as 
calculated by the FFS method with the FP times of the outer ones obtained from 
direct simulations. Furthermore we have proposed a numerical route to construct 
the arm relaxation correlation functions from the FFS simulation data saved on 
discrete interfaces. This method is essentially a summation of weighted exponential 
relaxation functions with characteristic times determined by the mean FP times of 
different slip-links along the arm. The so-constructed arm end-to-end vector cor- 
relation functions, ®(t), show reasonably good agreement with those obtained in 
standard slip-spring simulations, while larger quantitative discrepancies are found 
for the stress relaxation functions G(t) probably due to the biased selection of local 
segment conformations during the FFS simulations. 

We have also attempted to extend the FFS method to systems with constraint 
release, namely to entangled star polymer melts. The key change from the non- 
CR case is to define the reaction coordinate using the innermost surviving original 
slip-link. Again good agreement is found between the FFS simulation results on the 
terminal arm relaxation time with those obtained in standard slip-spring model sim- 
ulations. Therefore the combined FFS and slip-spring simulation method provides 
an efficient tool for studying the dynamics of highly entangled branched polymers 
which are generally inaccessible to direct simulation methods but highly desired for 


the development of quantitative theories on entangled branched polymers. 


Chapter 4 


Relaxation of Branched Polymers: A 
Combinational Study by Molecular 
Dynamics and Slip-Spring Model 


4.1 Overview 


To describe the rheological behaviour of branched polymers and their general mix- 
tures, the tube theory has to incorporate different relaxation mechanisms, such 
as contour length fluctuation or arm-retraction [27], and constraint release which is 
modelled by either dynamic tube dilation or constraint release Rouse motion 
\22|. A number of approximations and assumptions have to be made for describing 
various experimental results. For example, for describing 3-arm asymmetric stars, it 
was assumed that the full retraction of the short arm allows the branch point to hop 
a fraction of the tube diameter, pa, where p is a factor smaller than 1 and a is either 
the original or dilated tube diameter under different assumptions. In order to fit 
the experimental data, the fitting parameter p? ranges from 1 to 1/60 for different 
asymmetric stars. For H-polymers, however, the range of p” is relatively narrow, 


roughly from 1/12 to 1/15 [B6]. Recently Baéova et al. performed large-scale 
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molecular dynamics simulations of entangled branched polymers and found that 
considering hopping in the dilated tube provides the most consistent set of hopping 
parameters in different architectures. However, whether the value of p? should be 
universal or system-dependent remains unknown, which implies that both the the- 
oretical model and the underlying assumptions should be examined starting from 
microscopic principle. 

On the other hand, different theoretical or numerical models at more fine-grained 
levels have been developed for describing entangled polymers. Representative exam- 
ples are the slip-link based models [59H67]. This class of models treats the entangle- 
ments as binary contacts between different chain segments, and thus can introduce 
finer details, such as the conformation of polymers in space, the specified locations 
of entanglements, and the spectrum of constraint release rates. Among these mod- 
els, the slip-spring model developed by Likhtman has been shown to describe 
the MD simulations and experimental results on linear systems reasonably well [81]. 
Most recently, the slip-spring model has been extended to study symmetric stars and 
star-linear blends by comparing the stress relaxation modulus G(t) with the experi- 
mental data [10i). Cao and Wang have also used slip-spring model and the MD 
simulations based on Kremer-Grest (KG) model to investigate the arm-retraction 
and constraint release effects on star polymers, and examined the mechanism pro- 
posed by Shanbhag et al. on explaining the release of the deepest entanglements 
on the star arms. 

In this chapter, we investigate the application of the slip-spring model on branched 
polymers of different architectures, starting from testing the consistency of its pre- 
diction power, such as the universality of model parameters for both linear chains 
and branched polymers. The model systems we studied include 3-arm symmetric 
stars, asymmetric-stars, and H-polymers. For the later two architectures, current 
slip-spring model might fail because some mechanisms are missing, e.g., the slip- 


springs on the cross-bars of H-polymers cannot be released by arm retraction. In 
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order to find the additional mechanisms, microscopic understanding from molecular 
dynamics simulations is required. 

This chapter is arranged as follows. In Sec. we will introduce the advanced 
techniques for performing highly efficient MD simulations and data analysis of mildly 
entangled branched polymers represented by the fully flexible Kremer-Grest bead- 
spring model. The slip-spring model has two versions which are different on the way 
to handle slip-link motion. In the original version, the slip-links diffuse continuously 
along the chain backbone, following the standard Brownian dynamics (BD) [66]. In 
a recently updated version (see Sec. (1.6.2), the slip-links move discretely by hopping 
between neighbouring monomers as governed by a Monte-Carlo algorithm in order to 
achieve higher efficiency [83]. The optimization of the model parameters, including 
the frequency to perform MC algorithm f**, the coarse-graining parameter No, and 
the time-scale mapping factor to, and the test of their consistency between linear 
and branched polymers are given in Sec. In Sec. |4.4, we present the MD 
simulation results for symmetric stars, asymmetric stars, and H-polymers, whose 
observables, such as the end-to-end vector relaxation function ®(t) and the monomer 
mean-square-displacement g;(t) are compared with the predictions of the tube-based 
theory. Following that, we compare the simulation results obtained by running the 
standard slip-spring simulations with the MD data on the same systems and show 
the significant discrepancy between them, especially for asymmetric stars and H- 
polymers. Such discrepancy can be attributed to the missing mechanism mentioned 
before. To cope with it, a parameter-free algorithm allowing the slip-links to cross 
the branch-point will be added into the slip-spring model, whereby a remarkable 
improvement can be achieved on the agreement with MD results. The conclusions 


will be given in the last section. 
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4.2 MD Simulation Method 


The KG bead-spring model (see Sec. is the most widely used generic MD 
model for entangled polymers, in which the beads representing the monomeric units 
interact with each other via the purely repulsive Lennard-Jones (L-J) potential. 
Combined with the bonding potential modelled by FENE, the excluded volume 
interactions can effectively prevent the chains from crossing. Chain stiffness can be 
introduced into the KG model via a three-bead bending potential, whereby more 
entanglements can be implemented with same chain length. For distinction, the 
chain model with and without bending potential are, respectively, called fully flexible 
and semi-flexible KG models. The semi-flexible KG model is relatively cheaper on 
the computational cost, but the simulation results obtained using fully flexible KG 
chain model have been shown to have better agreement with experimental data; 
thus the fully flexible KG model is still the first choice as long as the computational 
cost is affordable |82} 115}. 

MD simulations of the entangled star polymers are extremely expensive, because 
the terminal relaxation time 74 grows exponentially with the number of entangle- 
ments per arm. In order to use the fully flexible KG model, we employed a high- 
performance GPU package called HOOMD [117 [118], which allows us to reach the 
terminal relaxation time of the symmetric stars with arm-length up to M = 384 
(for convenience we use M to represent the number of monomers in KG model 
and use N to represent the number of beads in slip-spring model in this chapter). 
Considering the average entanglement segment length is about 50 to 65 [82], 
the corresponding number of entanglements per arm is around 6 to 7. The sim- 
ulations are performed in the NVT ensemble with a periodic boundary conditions 
applied to all three dimensions. The volume of the central cubic box is obtained by 


V = Mr/p, where Miot is the total number of beads and p = 0.8507? is the density 
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of the beads. The equation of motion of beads is given by 


where r; is the coordinate of the i-th bead, m is the mass of beads, [ = 0.5(mkpT)/?/o 
is the friction coefficient, W is the Gaussian white noise satisfying (W;(t) - W,(t’)) = 
6,;0(t—t’)6kg ITI, and I denotes the three-dimensional unit matrix. The simulation 
time-step is At = 0.0127,, while 7,3 is the Lennard-Jones time. 

The equilibration of the initial system is carried out by a home-made code called 
generic polymer simulator (GPS). In order to achieve a faster equilibration, a 
soft potential [102] is employed to perform the relaxation. The potential functions 


for bonded and non-bonded interactions are formulated by 


where k, = 20€, ro = 1.2220, r, = 1.60, and up = 2.2€. The soft potential allows the 
chains to cross each other, but can preserve the static properties, such as the chain 
conformations, close to the KG model. 

The dynamic observables are obtained using the data analysis tools in the GPS 
code, which can efficiently calculate the time correlation functions on-the-fly via an 
algorithm called “correlator” [119]. During simulations, GPS can work as an exter- 
nal module of HOOMD. Specifically, HOOMD generates the trajectory coordinates 
of the beads, which are read and analyzed by GPS at high frequency; then, GPS 
update the data of measured observables, and stores the trajectory files at a low 


frequency for further analysis. 
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4.3. Calibration of the Slip-Spring Model Parame- 


ters 


4.3.1 Basic Parameters 


In Sec. we have introduced the slip-spring model that takes Monte-Carlo 
algorithm to govern the diffusion of slip-links. Slip-spring model requires three basic 
parameters: the average number of beads between neighbouring slip-links N°, the 
number of the beads on the virtual spring N°, and the frequency to perform the MC 
moves f*°*. In Ref. [66], Likhtman compared a variety of the combinations of NS* 
and N°* in the slip-spring model of the Brownian dynamics version, showing that 
their combinations determine the plateau modulus G, in stress relaxation. Thus the 
plateau regimes can be superimposed on each other by adjusting them in pair. In a 
standard setting, N°° and N5* are 4 and 0.5 respectively. Once N®° and N®* are 
decided, the entanglement segment length is determined, thus the coarse-graining 
level is also determined. Other pairs of N° and NS* could be employed to change 
the level of coarse-graining. For example, NS° = 8 and NS*° = 1 will reduce the level 
of coarse-graining by half, thus one must double the number of beads to preserve 
the number of entanglements. The finer slip-spring model shows finer resolution in 
early regimes with additional computational cost. 

The frequency of MC attempts per time-step f** is related to the friction co- 
efficient of the slip-links in the original BD version. & is supposed to have little 
effect on rheological properties when it is much smaller than the friction coefficient 
of beads &). In the standard setting of the BD version, €, is set to be 0.1€9 since 
technically it cannot be 0. Similarly, f° in the MC version must be large enough to 
ensure that the slip-link can efficiently find out the local minimum potential within 


each time-step. In order to find the optimal value of f°°, we compared the viscosity 
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Figure 4.1: The viscosity 7 obtained at different frequency f** in the slip-spring 
model for linear chain systems with NS° = 4 and N®* = 0.5. (a) The logarithm plot 
with f°* ranging from 0.1 to 50. (b) The linear plot zooming into the range of f** 
from 0.5 to 10. 


n obtained at different f°, where 7 is given by 


n= if G(t)dt. 


As shown in Fig. [4.1[a), we choose three different chain lengths, N = 25, 38 and 51, 
to test their viscosities at the frequencies f°* ranging over 2 decades, namely from 
0.1 to 50, with the time-step At = 0.0570 and 7 is the slip-spring unit time [66]. In 
Fig. [4.1[a), n decreases fast when f** is smaller than 1. Afterwards, it decays much 
slower with increasing f°°. If we zoom into the range from 0.5 to 10 and plot the 
data in the linear scale of {°° (see Fig. [4.1{b)), the decrement of the viscosities at all 
chain lengths has an obvious “plateau-like” region after {°° = 4. In fact, 7 decreases 
by 30% when f*° increases from 5 to 50. Similar phenomenon happens to the BD 
version, in which 7 drops by 30% when the friction of the slip-link decrease from 


0.1€ to 0.01€5 [66]. We conjecture such monotonic decease of 7 with increasing f** 
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or decreasing €* as a consequence that in finite time-step the slip-links affects the 
mobility of the monomers connected by them, and thus influence the dynamics of 
the whole chain. Therefore, this effect can be eliminated by adjusting the horizontal 


shifting factor or the time mapping factor to. 


9imia(t)/ re 
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Figure 4.2: The horizontally shifted middle monomer mean-square displacements of 
linear chains at different MC frequency f**. The chain lengths are N = 25, 38 and 
51 respectively. The shifting factors are adjusted to make the curves superimposed 


on each other at the chain length N = 25. 


For further verification, we compare the middle monomer mean-square displace- 
ments gi mia(t) of the chains at different f SS whose results are presented in Fig. |4.2 
The curves of gi mia(t) have been normalized by the Rouse power law t!/2 in order 
to bring the curves into one decade along vertical axis. Four frequencies are chosen, 
namely f°> = 0.1, 1, 5 and 50. According to Fig. it is expected that when 
f° > 4 the slip-spring model should exhibit the same dynamics in late regimes 
after adjusting the time mapping factor tg. For comparison, we shifted the curves 
of gimia(t) horizontally, making them superimposed at N = 25 and then use the 


same shifting factors for other chains lengths. As expected, the gi mia(t) curves for 
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f°* = 5 and 50 are superimposed in late regimes, while the curves for f°° = 0.5 
cannot overlap with those of other frequencies at the chain length N = 38 and 51. 
Therefore, we can define f° = 5 as a standard parameter of the MC version when 
At = 0.057. With smaller At, f° = 5 should be reduced proportionally, e.g., 
f°5 = 1 when At = 0.017. 


4.3.2 Mapping Parameters 


The parameters used to map slip-spring results to KG model data include the length- 
sale and time-scale mapping factors. On length scale, one must determine how many 
monomers in KG model are represented by one bead in the slip-spring model. This 
mapping number No has an unique value for a certain pair of NS° and N5°. On 
time scale, the exact solution of the mapping factor tp due to coarse-graining remains 
unknown. But its value must be a constant for a given model, thus can be easily 
found by fitting the dynamic observables, such as the end-to-end relaxation function 
®(t) or the stress relaxation function G(t). In Ref. [82], Wang et al. explored the 
parameter sets. When NS* = 4 and N°* = 0.5, they found the mapping parameters, 
to = 3370 and Nop = 9.74 for fully flexible KG linear chain model. On length scale, 
the mapping factor is non-arbitrary, but determined by No and C.,, where Cy, is 
the characteristic ratio of the chain. For example, the vertical shifting factor for the 
end-to-end vector relaxation function ®(t) is the product of C,, and No. In fully 
flexible KG model, C,, is around 1.82 [5I]. 

It is expected that the length-scale mapping factor is independent of f°°, but 
the time mapping factor ty can be affected by f°° according to our previous discus- 
sion. To get the value of tg, we performed a series of MD simulations using fully 
flexible KG model of linear chains, whose lengths are M = 256, 320, 384, 448 and 
512, respectively. With No = 10, the chain lengths used in the slip-spring model 
simulations are N = 26, 32, 38, 45 and 51 respectively. As shown in Fig. [4.3{a), 
the data of the slip-spring model on ®(t) agree well with the MD results at all chain 
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Figure 4.3: Mapping of the slip-spring model results (lines) obtained with fS° = 5 
and At = 0.0579 to the data of fully flexible KG model (symbols) for linear chains on 
the end-to-end relaxations ®(t) and the middle monomer mean-square displacements 


1,mia(t). 


lengths when using tp = 3400. This value is very close to the previous work [82]. 
A further check is carried out on the middle monomer mean-square-displacements, 
gimia(t), as shown in Fig. [4.3(b). This is a more strict examination, because the 
single bead diffusion is very sensitive to the slip-spring parameters, especially after 
T,. Again, g1mia(t) is divided by t'/? to bring the curves into one decade for better 
comparison. The standard parameter setting fits the simulation data very well in 
the middle and late regimes. The fitting in early regime, t < T., could be improved 
by finer-graining. In Refs. and [82], the finer-grained slip-spring model simu- 
lations were performed with N°S = 8 and NSS = 1, which gives better resolution in 


early regimes. 
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4.4 Relaxation of the Branched Polymers 


4.4.1 Simulation Systems 


Fig. |4.4] presents the schematic plot of the branched polymer architectures inves- 
tigated in this work, including symmetric stars, asymmetric stars and H-polymers. 
For convenience, several subscripts are added on “N” and “M” to denote different 
architectures, i.e., “sym” for symmetric star, “asy” for asymmetric star, and “h” for 
H-polymer. The superscripts are added on “N” and “M” to denote the types of 
subchains, i.e. “I” and “s” represent long and short arm in an asymmetric star, 
while “t” and “a” respectively represent the cross-bar and the arm in a H-polymer, 


respectively. 


Masy Me 
M. sym 


Ming 
(a) (b) (c) 


Figure 4.4: Sketches of branched polymer architectures: (a) Symmetric star, (b) 


Asymmetric star, (c) H-polymer. 


The MD simulation systems are listed in Table Five symmetric stars are 
investigated, whose arm lengths vary from 64 to 384. For asymmetric stars, two 
long-arm lengths are chosen, i.e., Mie = 256 and 384. We only choose one cross-bar 
length M;} = 256 in H-polymers due to extraordinarily long relaxation time. In 
all cases, there are 100 molecules in the central simulation box. Thus the biggest 
system has 115,300 particles. The number of molecules is large enough to ensure 
the end-to-end distance or radius of gyration is less than 2/3 of the cubic box size, 


whereby the finite size effect is negligible. With a time-step At = 0.0127,), the total 
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Symmetric Star Asymmetric Star H-Polymer 
Mgym  Neym || Masy Mésy | Nasy Nasy |] Mn Ma | Na Net 
64 6 256 = 64 25 6 256 64 | 25 6 
128 13 256 128 | 25 13 | 256 128] 25 13 
256 25 384 364 38 6 

320 a2 384 128 | 38 13 

384 38 384 256 | 38 20 


Table 4.1: Molecular structure parameters used in the KG model (/), and the 


corresponding slip-spring model (NV). 


simulation time reaches 5 x 10’7,3, covering a time range of 9 decades. 

The simulation systems of slip-spring model are also listed in Table In 
the slip-spring model, each bead corresponds to 10 monomers in the KG model. 
Accordingly, each branch-point bead in the 3-arm stars or H-polymers represents the 
corresponding branch-point monomer in the KG model plus 3 monomers from each 
subchain connected to it. In all slip-spring model simulations, the total number of 
molecules in the ensemble is set to be 30, which is sufficient to obtain good statistics 
of the measured observables. Each simulation runs roughly about 2 hours on a single 


CPU to achieve more than 10 terminal relaxation times. 


4.4.2 MD Simulation Results 


In this subsection, we will focus on the monomer mean-square displacements of dif- 
ferent branched polymers obtained from MD simulations. Those observables provide 
rich microscopic information in microscopic dynamics, and are usually compared 
with the relaxation regimes predicted by the tube theory [17]. 

Fig. [4.5] presents the middle monomer mean-square displacements gj mia(t) of the 
symmetric stars. In Fig. [4.5{a) and (b), the gi. mia(t) data (open symbols) have been 
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(a) (b) 
Figure 4.5: MD results on the mean-square displacements of the middle monomers 
of the arms g; mia(t) (open symbols) and the branch points 91 prancn(t) (solid symbols) 


for symmetric stars. 


divided by t!/? and t'/4 to reveal different regimes, respectively. At time scales t < Te, 
the monomers are not aware of the topological constraints imposed by surrounding 
chains, and thus follow the standard Rouse motion, gi mia(t) ~ ¢/?. As a result, 
the first plateau-like regime is found before 7, (dashed line) in Fig. [4.5{a). At time 
scales T. < t < Tg, the monomer diffusion follows “constraint Rouse”, where the star 
arm experiences 1D Rouse relaxation in the confining tube, leading to a power law of 
g(t) ~ t/4. Thus, in Fig. [4.5{b), a plateau regime can be found in gi mia(t) curves 
after T., whose width increases with the growing arm length due to the power law 
of TR ~ Me. For Mgym = 64, the barely observed plateau implies that Tp of the 
arm is roughly equal to 7 and thus the entanglement segment length is about 64, 
which is consistent with the measured M, & 50 — 65 by mean-square displacement 
data in Ref. [82]. In tube theory, a scaling law of gi(t) ~ t'/* is predicted at time 
scales TR < t < Tg. However, this region is hard to observe in linear chain melts 


|82|, because 7g and 7g are proportional to the square and cubic of the chain length 
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respectively, which requires very long chains to distinguish Tp and 7, in logarithm 
scale. The g, ~ t'/* power law is predicted for reptation (1D random walk in the 
tube) of linear chains. For stars, the relaxation proceeds by arm retraction, the 
behaviour could be somewhat different. But the separation between Tp and Tg due 
to the exponentially slow relaxation leads to a regime where g;(t) grows slowly after 
Tr. At time scales t > 74, the monomers follow free diffusion, gi mia(t) ~ t. It is 
important to note that the semi-flexible KG model leads to different power laws in 
two regimes, t < T. and T. < t < TR, where g;(t) are proportional to t°° and t®? 
respectively [82]. Therefore, the fully flexible model agrees better with the tube 
model. 

Another interesting observable is the branch point mean-square displacement 
Ji,branch(t). The gi branch(t) data for symmetric stars are shown by the solid symbols 
in Fig. |4.5| Due to the connections with more than 2 monomers, the Rouse motion 
of the branch point is different from other monomers. At time scales T. < t < Tr, 
Ji,branch(t) roughly follows t'/5 scaling for Msym > 256 rather than t'/4, due to the 
cage effect. In Fig. [4-5[a), the minimum of the gi branch(t) curve corresponds to 
the terminal relaxation time Tg, after which the branch point follows free diffusion, 
G1,branch(t) ~ t. 

For asymmetric stars, the 9) branch(t) exhibits a significant speeding-up due to the 
shortening of one arm. Fig. [4.6{a) presents a comparison between the 91 branch(t) 
data of the symmetric and asymmetric stars. The arm lengths of the symmetric 
stars are 256 and 384. Shortening one arm of the symmetric stars into 64, 128, or 
256 for Nsym = 384, the mobility of the branch point increases significantly at time 
scales T, < t < Tg. The earlier terminal relaxation brought by shortening one arm is 
reflected in the arm end-to-end relaxation ®(t). As shown in Fig. [4.6{b), we present 
®(t) of long arms of the asymmetric stars, which are compared with ®(t) of the 
symmetric stars. The linear chains with the lengths of 512 and 768 are treated as 


2-arm symmetric stars, whose arm end-to-end relaxations ®,,;q(t) are also plotted in 
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Figure 4.6: (a) MD results on branch point mean-square displacement g1 branch (t) of 
the symmetric and asymmetric stars, and middle monomer mean-square displace- 
ment gj mia(t) of the linear chains. (b) End-to-end relaxation ®(t) of the arms of 
the symmetric stars and the longer arms of the asymmetric stars, and middle-to-end 


relaxation ®,,iq(t) of the linear chains. The same symbols are used in both figures. 


Fig. [4-6{b). The terminal arm end-to-end relaxation time Tg of the symmetric star 
with the arm length Mj, = 256 is around 3.8 x 10°r,3. For the asymmetric star 
with Mey = 256 and M&.. = 64, Tq of the longer arms is around 2.5 x 10°73, which 
is almost equal to Tg of the linear chains with a length of 512. In this case, the short 
arms seem to have little effect on the terminal relaxation time when attaching them 
to the middle monomer of linear chains, but leads to a slower Rouse motion due to 
extra connection. 

For H-polymers, we present the mean-square displacements of the branch points 


Jibranch(t) and the middle monomers of the cross-bars gj ,,iq(t), as shown in Fig. 
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Figure 4.7: MD results on mean-square displacement of the branch points 91 prancn(t) 
and the middle monomers gj ,iq(t) in the H-polymers together with that of middle 


monomers of linear chains. 


4.7| The normalized 91 branch(t) reaches the minimum earlier than Ji mia(t) of cross- 
bar by more than a decade, implying that the arm relaxation happens much earlier 
than the cross-bar relaxation. According to the “hierarchical” hypothesis [39] {120}, 
after the relaxation time of the branch arms, the cross-bar do reptation with higher 
effective frictions at the two ends arisen from the relaxed arms. Therefore, the shape 
of J} mia(t) of cross-bar should be similar to that of linear chains. In Fig. one 
can find that gi mia(t) of the linear chain with a length of 512 is almost the same as 
Ji mia(t) of the H-polymer with Mj = 256 and Mf = 64. It could be a coincidence 
since the molecular weight of them are equal, but can also imply that the short arm 


length M;’ = 64 has weak entanglement. 


4.4.3 Slip-Spring Model for Branched Polymers 


In this subsection, we first present the simulation results of the current slip-spring 


model with MC moves on the branched polymers discussed above. Fig. [4.8{a) 
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presents the results of symmetric stars, including the arm end-to-end relaxation ®(t), 
the branch point mean-square displacement gi pbranch(t) and the middle monomer 
mean-square displacement of arms gi mia(t). The slip-spring model results have 
been shifted to map the MD simulation data with the mapping parameters exactly 
the same as those used in linear chains, i.e., No = 10 and ty = 3400. It is found 
that the agreement is good for all observables. Therefore, current slip-spring model, 
which has only CR and arm-retraction as the relaxation mechanisms, is sufficient to 
describe the slightly and mildly entangled symmetric stars. We can expect that it 
also works for the well-entangled symmetric stars. This is reasonable, because for 
such polymers, the terminal relaxation times of the systems are the same as the arm 
retraction times. 

The simulation results of the asymmetric stars are given in Fig. [4.8{(b), including 
the arm end-to-end relaxation ®(t), the middle monomer mean-square displacement 
of long arms G1 mia () and the branch point mean-square displacement 9} branch(t). In 
these plots, significant discrepancies can be found in the 9) branch and 1 mia(t) data 
for most asymmetric stars except the one with Macy = 384 and Mj. = 256. We 
presume this is due to the lower asymmetricity, i.e., the short arm length is relatively 
closer to the long arm length than other asymmetric stars. The discrepancies of 
other asymmetric stars occur at time scales close to Tq of long arms, where 9} »niq(t) 
and gi,branch(t) curves predicted by the slip-spring model are both lower than the 
MD results. The significantly underestimated mobility of branch points and middle 
monomers indicates that CR and arm-retraction are insufficient to describe the 
relaxation of asymmetric stars, especially for those with larger asymmetricity. On 
®(t), however, the slip-spring results agree well with the MD results for both short 
arms and long arms, which might be because the shape of ®(t) is not strongly 
affected by the branch point motion. 

In H-polymers, the current slip-spring model also does not work. As shown in 


Fig. |4.9| @(t) of the cross-bar reaches a plateau after the Rouse relaxation, because 
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Figure 4.8: Simulation results of the previous slip-spring model and the KG model 
on the end-to-end relaxation ®(t), the branch point mean-square displacement 
Jipranch(t), and the middle monomer mean-square displacement 9) mia(t) for (a) 
the symmetric stars, and (b) asymmetric stars. In asymmetric stars, we only plot 
9imia(t) for long arms. The symbols and lines are the results of the KG model and 
the slip-spring model, respectively. In bottom plot of (b), the solid symbol and lines 


are for the long arms, while the open symbols and dashed lines are for the short 


arms. 
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Figure 4.9: The simulation results of the slip-spring (lines) and MD using KG model 
(symbols) on ®(t) of the H-polymers. 


the slip-links on the cross-bar are blocked at the two branch points. The results in 
Fig. [4.8{b) and imply that the failure of slip-spring model on asymmetric stars 
and H-polymer may originate from the same problem: the entanglements should be 
allowed to cross the branch point in some conditions; or in other words, the branch 
point can “hop” along the confining tube of the two long arms of an asymmetric 
star, or the confining tube of the cross-bar and one arm in a H-polymer. 

In Ref. [63], Shanbhag and Larson proposed a slip-link model for branched 
polymers, in which they assume that the branch point hops when the slip-links on 
the short arm are all removed, corresponding to a complete arm retraction event. 
The advantage of this assumption is that it requires no additional parameters. In this 
work, we can introduce a parameter-free slip-link “hopping” mechanism in current 
slip-spring model, i.e., the slip-link can cross the branch point and hop onto another 
arm when the 3rd arm is fully relaxed. This hopping step is taken using the standard 
MC step of slip-spring model, thus requiring no additional parameters. It must be 
noted that the slip-link hopping in H-polymers only occurs between cross-bar and 


one arm. 
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With this slip-link “hopping” mechanism, we plot the results of the slip-spring 
model for symmetric stars in Fig. {4.10} The results on arm end-to-end relaxations 
@(t) and mean-square displacements for both branch points g1 pranch(t) and middle 
monomers of arms gj mia(t) show little difference from Fig. [4.8[a). In Fig. [4.10{d), 
we also present the average waiting time of one “hopping” event for each subchain, 
Thop, Which is proportional to the inverse of the hopping rate on each branch point 
Ties 


Thop = ats, (4.2) 


Thop 

where g is branch point functionality (i.e. the number of subchains connected to one 
branch point). In one time-step, a slip-link can cross a branch point many time when 
f* > 1, we consider the frequent crossings due to the algorithm of slip-link diffusion 
should only be counted once for each time-step if the hopping succeeds. As shown in 
Fig. [4.10(4), Thop 18 equal to Tg of the arm end-to-end relaxation function when the 
arm length Msym is smaller than 128, after which Top is even larger than 7g. Thus, in 
symmetric stars, contribution of slip-link hopping to relaxation is negligible, which 
also explains why the original slip-link model can well describe symmetric stars. 

For asymmetric stars, the relaxation times of the short arms 73 and the long 
arms T} are well separated, thus the relaxation can be accelerated by the “hopping” 
mechanism. Fig. [4.11{a) shows the mean-square displacements of branch points 
Ji,branch(t) and middle monomers for both long arms gj q(t) and short arms g? jy.q(t) 
in the slip-spring model simulations, which are found to be in good agreement with 
the MD results of KG model. Comparing to Fig. [4-8{(b), the slip-spring model 
with slip-link hopping mechanism not only predicts the separation of Gh mialt) of the 
asymmetric stars with different short-arm lengths, but also well predicts the higher 
mobility of branch points in the terminal regime. In Fig. [4-11{b), the arm end-to- 
end relaxation ®(t) is similar to Fig. [4.8{b), showing that the hopping mechanism 
does not affect the prediction of &(t). In Fig. [4-11{c), Thop Of the asymmetric stars as 


calculated using Eq. [4.2]are compared with the end-to-end terminal relaxation times 
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Figure 4.10: Simulation results of the slip-spring model with slip-link “hopping” 
mechanism at the branch point in comparison with the MD data using the KG model 
for the symmetric stars: (a) the branch point mean-square displacements 9) pranch(t), 
(b) the middle monomer mean-square displacements gi mia(t), (¢) the arm end-to- 


end vector relaxations ®(t), (d) the average waiting times of one “hopping” event 


for each subchain 7Tpop and the end-to-end terminal relaxation times Tg. 
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Figure 4.11: Simulation results of the slip-spring model with slip-link “hopping” 
mechanism at the branch point in comparison with the MD data using the KG model 
for the asymmetric stars: (a) the mean-square displacements of the branch points 
Ji,branch(t), the middle monomers of long arms gj ,iq(t), and the middle monomers 
of short arms gj ,ia(t), (b) the arm end-to-end relaxations (¢) (solid symbols and 
lines represent long arms, open symbols and dashed lines represent short arms), (c) 
the average waiting time of one “hopping” event for each subchain top, and the 


end-to-end terminal relaxation times of short arms 7$ and the long arms 7}. 
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Figure 4.12: Simulation results of the slip-spring model with slip-link “hopping” 
mechanism at the branch point in comparison with the MD data using the KG model 
for the H-polymers: (a) the branch point mean-square displacements 91 branch(t), (b) 
the mean-square displacements of the middle monomers of the arms g?',,;q(¢) and 
the cross-bars gj ia(t), (c) the end-to-end relaxations ®(t) (solid symbols and lines 
represent cross-bars, open symbols and dashed lines represent arms), (d) the average 
waiting time of one “hopping” event Thop for each subchain, and the end-to-end 


terminal relaxation times of the arms 7 and the cross-bars T}. 
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of long arms 7} and short arms 7§. With identical short-arm length, 7),.) for different 
long-arm lengths are equal, indicating that the relaxation of short arms dominates 
the “hopping” probability. For most asymmetric stars, Thop is smaller than 75 apart 
from the asymmetric star with Mi = 384 and M2,, = 256, whose Thop is equal to 
Ti. Thus we can affirm that the hopping mechanism is trivial in the systems with 
low asymmetricity, which also explains why the slip-spring model without “hopping” 
mechanism can well describe this asymmetric star (see Fig. [4.8{b)). 

In H-polymers, the cross-bar relaxation is dominated by the diffusion of branch 
points, thus a significant improvement of the slip-spring model for H-polymers can 
be achieved by incorporating such slip-link “hopping” mechanism, as shown in Fig. 
4.12} For all observables considered, the agreement with the KG model is reasonably 
good. In Fig. [4.12(a), the gi branch(t) data for M? = 64 is higher than the MD results 
by roughly 20% at time scales 7. < t < Tg. Similar disagreement can be found in the 
middle monomer mean-square displacements of the arms, as shown in Fig. [4.12{b). 
The overestimated branch point mobility for M@? = 64 is also reflected in the ®(t) 
plots in Fig. [4.12{c), where the relaxation of the cross-bar is slightly faster than 
the MD result. This might result from the fast creation and deletion of the single 
slip-link at the free ends of the short arms which on average have less than two 
slip-springs. For H-polymer with longer arm length M;* = 128, the agreements on 
these observables become much better. In Fig. [4.12{d), Thop for both H-polymers 
are much shorter than the end-to-end terminal relaxation times of the arms 7”. 
Thop for My = 64 and M;) = 128 are roughly half of 7.) for asymmetric stars with 
Misy = 64 and M3., = 128, because there are two short arms connected to each 


branch point in H-polymers. 
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4.5 Conclusions 


In this chapter, we presented a detailed description of the slip-spring model using 
Monte-Carlo algorithm to govern the motion of the slip-links. This version of the 
slip-spring model enables us to employ a larger time-step than the original Brownian 
dynamics version without changing the predicted dynamic properties. With a careful 
selection of the parameters, we obtained the simulation results of the slip-spring 
model on monodisperse linear chains, which are found to be in good agreement with 
those of the molecular dynamics simulations using the fully flexible Kremer-Grest 
model. In a standard parameter setting, i.e., NSS = 4 and NS* = 0.5, the mapping 
parameters for the fully flexible KG model are found to be ty) = 3400 for time scale 
and No = 10 for length scale, which are believed to be universal for different polymer 
architectures. 

After the careful calibration, the slip-spring model was extended to branched 
polymers. In order to examine the slip-spring model results, we performed exten- 
sive molecular dynamics simulations using the flexible KG model for a variety of 
branched polymers, including the 3-arm symmetric and asymmetric stars, and the 
H-polymers. The slip-spring model results agree well with MD simulations on sym- 
metric stars, but fail on both asymmetric stars and H-polymers. We consider this 
problem originating from the absence of some relaxation mechanisms in the model. 
One possible mechanism is that the slip-links can hop over a branch point between 
two subchains connected to it when the third arm has no slip-links on it. The 
modified slip-spring model with the “hopping” mechanism provides results in good 
agreement with MD data without requiring extra parameters. Such “hopping” mech- 
anism is found to be important for asymmetric stars with significantly different arm 


lengths and essential for H-polymers. 


Chapter 5 


Conclusions 


The study of the dynamics of entangled branched polymers is of both fundamental 
and practical importance. In the thesis, we focus on investigating the relaxation 
behaviours of entangled branched polymer with simple architectures. 

Due to the steeply growing quadratic potential, arm retraction is a typical mul- 
tidimensional first-passage problem, whose exact solution remains a open question. 
In the Milner-McLeish theory [30], this problem is simplified by treating the whole 
chain as one bead attached to the branch point via a harmonic spring. This so- 
lution, however, overestimated the relaxation time by neglecting the contributions 
of multiple Rouse modes. Cao et al. presented an analytical solution to the 
multi-dimensional first-passage problem, which predicts an arm relaxation time in 
the absence of constraint release 2/N times smaller than that given by the Milner- 
McLeish theory, where N is the arm length. In order to examine the two theoretical 
models, we employ advanced numerical methods for studying first-passage problems, 
such as the forward flux sampling and weighted ensemble methods. These methods 
were first implemented to the simplest multi-dimensional first-passage problem: the 
escaping time of a Brownian particle in a 2D potential well. Due to its remarkable 
performance in all aspects, the FFS method was then chosen to study the exten- 


sion problem of a 1D Rouse chain model, which is analogous to the arm-retraction 
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problem. We found that the first-passage time is getting shorter if the Rouse chain 
is represented by more beads, showing good agreement with the prediction of the 
asymptotic solution of Cao el al. [54]. 

Then the FFS method was implemented to solve the arm-retraction problem in 
the slip-spring model, which is a coarse-grained bead-spring model for entangled 
polymers. With a controllable precision, this method allows direct comparison be- 
tween the slip-spring model and the tube theory for well-entangled star polymers 
with up to 16 entanglements per arm. Moreover, a study is conducted on the ex- 
traction of experimentally measurable observables from FFS simulations, such as 
the end-to-end vector and stress relaxation functions. We believe this work will not 
only expand the application of FFS method to polymer dynamics by reproducing 
full dynamic spectrum rather than just the first-passage time, but also to many 
other scientific areas. 

After the remarkable success on linear melt systems, the slip-spring model has 
been extended to the study of the branched polymers [34] [LOI]. However, due 
to the absence of certain mechanisms, current slip-spring model cannot describe 
the relaxation behaviours of some architectures, such as asymmetric stars and H- 
polymers. To cope with it, we conducted a series of MD simulations for branched 
polymers with different architectures, including 3-arm symmetric and asymmetric 
stars, and H-polymers. With the fully flexible Kremer-Grest chain model, these 
simulations achieved the terminal relaxation time of the mildly entangled arms with 
up to 7 entanglements. The slip-spring model, whose parameters have been carefully 
calibrated according to the MD results of linear chains, was implemented to predict 
the relaxation behaviours of these branched polymers. Comparing the MD and 
slip-spring model simulation results, we found a significant discrepancy due to the 
missing mechanism aforementioned. In the tube theory for asymmetric stars, it is 
assumed that the branch point can hop a fraction of tube diameter when the short 


arm is fully relaxed. For the slip-spring model, we proposed a slip-link “hopping” 
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mechanism, which allows the slip-links to cross the branch point when the third 
arm is unentangled. Using this mechanism, the slip-spring model is shown to have 
a good agreement with the MD results. 

In conclusion, many problem of entangled branched polymers dynamics remains 
open. Based on the achievements made in this thesis, future studies can be carried 
out in many aspects. One of the main challenges of branched polymer dynamics 
is the lack of clear microscopic picture of entanglements, which leads to numer- 
ous assumptions and fitting parameters in the current tube model. One solution 
is mapping a more fine-grained model onto a more coarse-grained one. This could 
effectively avoid the dependence on the exact definition of entanglements, while the 
missing mechanisms can still be revealed. In fact, this approach has been imple- 
mented in the fourth chapter to develop the slip-spring model for branched poly- 
mers by comparing its main results with MD simulations. In the second and third 
chapters, we have shown that the analytical and numerical solutions of the multi- 
dimensional first-passage time problem are effective to obtain the relaxation spectra 
of well entangled symmetric stars, which can be either modelled by the one dimen- 
sional Rouse model or the slip-spring model. Both models can be mapped onto the 
tube model, but further more their results can be compared with experiments and 
MD simulations, which could be useful to refine the tube theories, especially with 
constraint release. Another solution is to investigate the dynamic behaviours by 
analysing the evolution of entanglements in MD simulations. The powerful analysis 
tools, such as the “primitive path analysis” [51], and the newly developed “contact 
map analysis” and “tube axis” can pave the way to decipher the hidden 
mechanisms. For example, the constraint release events can be monitored by trac- 
ing the creation and destruction of entanglements in the middle of the chains, and 
the possible entanglement hopping mechanism of branched polymers can be detected 
by tracing the movements of entanglements close to the branch points. From these 


analysis, the distribution of the constraint release rates and the hopping distance 
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are available. Moreover, such an approach is very important to the investigation of 
entanglements in the nonlinear regime, since the solid evidences are most likely to 
come from MD simulations. In particular, one can figure out how does the num- 
ber of entanglement change under flow, and thus help to improve the theories and 


models for nonlinear dynamics. 
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